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ABSTRACT 


Weirs,  Vincent  Gregory.  Simulation  of  Ram- Accelerator  Flowfields  (Under  the 
direction  of  Dr.  Graham  V.  Candler) 

Ram-accelerators  are  devices  that  use  combustion  to  accelerate  projectiles  to  high 
velocities.  Shaped  projectiles  are  propelled  through  an  accelerator  tube  by  expansion  of 
combustion  products.  Ignition  of  the  premixed  fuel  and  oxidizer  is  shock-induced.  Tailor¬ 
ing  the  shape  of  the  projectile  results  in  improvements  in  performance.  In  this  research 
numerical  simulations  of  the  reacting  flowfield  surrounding  the  projectile  are  obtained  and 
analyzed. 

Computational  Fluid  Dynamics  (CFD)  methods  are  used  to  solve  the  extended  Navier- 
Stokes  equation  set.  Chemical  kinetics  are  modelled  by  a  quasi-global  reaction  mechanism 
including  10  species  and  12  reactions.  Flow  is  assumed  to  be  laminar.  The  governing 
equations  are  discretized  using  a  finite-volume  formulation.  Inviscid  fluxes  are  calculated 
using  a  Roe-TVD  scheme,  while  standard  central  differences  are  used  for  the  viscous 
terms.  A  modified  form  of  the  diagonal  implicit  method  is  used  to  advance  the  solution 
from  initial  conditions  to  the  steady  state.  Some  species  mass  conservation  equations  are 
replaced  by  elemental  mass  conservation  equations  to  improve  convergence.  As  a  whole, 
the  simulations  demonstrate  the  value  of  CFD  to  the  aerospace  engineering  field.  Physical 
models  are  evaluated,  numerical  methods  are  tested,  and  flow  phenomena  are  analyzed. 

In  the  course  of  the  research,  the  advantages  and  drawbacks  of  the  quasi-global  mech¬ 
anism  are  exhibited.  Compared  to  a  more  detailed  model  the  reaction  front  is  smeared,  and 
the  model  must  be  calibrated  with  established  data.  However,  the  quasi-global  mechanism 
offers  reasonable  accuracy  at  a  significantly  lower  computational  cost  than  more  detailed 
models. 

More  research  must  be  devoted  to  the  modified  diagonal  implicit  scheme.  For  some 
choices  of  replaced  species,  the  method  allows  negative  mass  fractions  to  develop  in  the 
evolution  of  the  simulations.  These  negative  fractions  can  persist  and  may  remain  in  the 
obviously  aphysical  steady  state  solution.  Careful  testing  was  done  to  ensure  the 


difficulties  were  avoided  in  the  results  presented  in  this  thesis.  However,  a  thorough 
inquiry  must  be  made  to  expose  the  underlying  problems  with  the  modified  diagonal 
implicit  scheme. 

The  ram-accelerator  simulations  assess  the  performance  of  the  double-cone  projectile 
geometry.  The  double-cone  forebody  was  proposed  to  delay  unstarts  to  higher  Mach  num¬ 
bers.  No  unstarts  were  observed  for  Mach  numbers  ranging  from  7  to  15.  The  simulations 
show  the  projectile  performs  as  designed.  A  large  separated  region  of  flow  is  located 
behind  the  projectile  shoulder.  It  is  the  result  of  a  shockwave/boundary  layer  interaction 
on  the  afterbody.  The  size  of  the  separation  is  affected  by  the  ignition  location.  The  sepa¬ 
ration  is  largest  when  combustion  and  consequent  energy  release  occurs  only  downstream 
of  the  region;  for  that  situation  the  expansion  of  combustion  product  enhances  the  shock¬ 
wave/boundary  layer  interaction.  Further  numerical  studies  are  necessary  to  assess  perfor¬ 
mance  for  different  fill  pressures  and  gas  mixtures.  More  detailed  physical  models  could 
also  be  included. 
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CHAPTER  I 


BACKGROUND  AND  MOTIVATION 


1.1  Introduction 

A  ram-accelerator  is  a  device  for  accelerating  projectiles  to  high  velocities.  It  relies  on 
chemical  energy  release  in  a  thermodynamic  cycle  similar  to  that  of  a  ramjet  engine.  Prac¬ 
tical  research  and  development  has  been  underway  for  less  than  a  decade.  While  the  con¬ 
cept  has  been  proven,  there  are  still  many  technical  issues  to  resolve.  The  operation  cycle 
can  break  down  for  several  reasons,  and  the  projectile  consequently  looses  thrust. 

Numerical  simulations  are  valuable  to  the  ram-accelerator  and  computational  fluid 
dynamics  (CFD)  communities.  To  the  former,  information  difficult  or  impossible  to  mea¬ 
sure  experimentally  is  made  available.  In  this  thesis,  the  conditions  which  lead  to 
‘unstarts’  are  sought.  CFD  benefits  from  the  validation  of  numerical  models;  this  will 
become  apparent  in  the  following  chapters. 

In  this  chapter  the  ram-accelerator  is  described;  the  theory  underlying  its  operation  and 
the  applications  of  such  a  device  are  presented.  Several  proposed  unstart  mechanisms  are 
described.  During  the  course  of  the  discussion  the  experimental  and  numerical  investiga¬ 
tions  to  date  are  reviewed.  Finally,  the  objectives  of  this  thesis  are  outlined. 

1.2  The  Ram-Accelerator  Concept 

The  ram-accelerator  concept  details  a  gasdynamic  method  of  accelerating  projectiles. 
The  components  consist  of  a  long,  constant-diameter  tube  and  a  shaped  projectile.  Initially 
the  accelerator  tube  is  filled  with  premixed  rations  of  fuel,  oxidizer,  and,  in  most  cases,  a 
diluent.  The  mixture  is  quiescent  and  pressurized  to  between  20  and  50  atm.  The  projectile 
at  its  widest  has  a  cross-sectional  area  less  than  that  of  the  tube.The  forebody  is  generally 
conical.  The  afterbody  usually  tapers  to  a  circular  base.  Fins  on  the  afterbody  align  the 
projectile  with  the  tube. 

The  projectile  enters  the  tube  at  a  supersonic  velocity,  generally  between  Mach  2  and 
Mach  4;  an  external  means  of  acceleration,  such  as  a  traditional  powder  gun  or  light  gas 
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gun,  is  used  as  a  preaccelerator.  Once  the  projectile  enters  the  accelerator  tube,  an  oblique 
shock  forms  off  the  nose.  It  reflects  off  the  tube  wall  and  the  projectile,  compressing  and 
heating  the  gas.  Eventually  these  shockwaves  raise  the  temperature  enough  to  ignite  the 
fuel,  and  chemical  energy  is  released  into  the  flow.  As  the  burned  gas  expands,  it  acceler¬ 
ates  the  projectile.  A  schematic  is  shown  in  Figure  1.1. 


Figure  1.1  Schematic  of  a  generic,  ram-accelerator  projectile  flowfield. 

The  accelerator  tube  is  initially  pressurized,  with  thin  diaphragms  sealing  the  ends. 
Usually  the  tube  is  composed  of  separate  sections,  which  can  be  individually  pressurized 
with  different  gas  mixtures.  By  staging  the  fill  pressures  and  mixtures,  the  acceleration  of 
the  projectile  can  by  controlled. 

Three  modes  of  ram-accelerator  operation  have  been  distinguished.^  They  are  based 
on  the  speed  of  the  projectile  relative  to  the  Chapman- Jouget  (C-J)  detonation  wave  speed. 
The  subdetonative  regime  extends  up  to  90%  of  the  C-J  velocity.  Combustion  is  subsonic 
(in  the  projectile  frame  of  reference)  and  occurs  primarily  in  the  wake.  One  or  several  pro¬ 
jectile  lengths  downstream,  the  flow  is  thermally  choked.  The  transdetonative  regime 
extends  from  90%  to  1 10%  of  the  detonation  wave  speed;  combustion  may  be  of  subsonic, 
supersonic,  or  mixed  character.  In  the  superdetonative  regime  the  projectile  velocity  is 
above  110%  of  the  C-J  speed,  and  generally  combustion  is  supersonic.  Note  that  these 
regimes  are  only  roughly  defined,  and  other  modes  have  been  proposed. 

1.  Hinkey,  Burnham,  and  Bruckner  (1992)  describe  these  operational  modes. 

2.  See,  for  example,  Hertzberg,  Bruckner,  and  Bogdanoff  (1988). 


Theoretically  projectiles  can  be  accelerated  to  velocities  as  high  as  10  km/sec.  So  far 
velocities  to  ~  2.7  km/sec  have  been  obtained  in  the  38  mm-diameter  tube  at  the  University 
of  Washington,^  which  has  the  oldest,  most  developed  ram-accelerator  research  program. 
The  majority  of  experimental  work  has  been  on  the  subdetonative  regime.  Since  the  dem¬ 
onstration  of  the  concept,  experimental  research  has  been  initiated  at  the  U.  S.  Army 
Research  Laboratory  (ARL),^  where  a  120  mm  accelerator  has  been  fired,  and  at  the 
French-German  Research  Institute  of  Saint-Louis  (ISL),^  where  a  90  mm  ram-accelerator 
is  being  developed. 

The  flowfield  within  the  ram-accelerator  tube  is  complex.  Besides  the  series  of  shock- 
waves,  there  are  contact  surfaces,  expansions,  boundary  layers,  flow  separations,  chemical 
induction  zones,  and  energy-release  fronts.  All  of  these  features  interact  with  each  other, 
further  complicating  the  analysis.  Numerical  simulations  must  accurately  resolve  these 
flow  phenomena  and  interactions.  Computational  flowfield  solutions  carried  out  at  ARL^ 
show  that,  in  general,  a  nonequilibrium  treatment  of  combustion  chemistry  is  required.  At 
Science  Applications  International  Corp.  (SAIC)^  intermolecular  force  effects  (real  gas 
effects)  have  been  shown  to  be  significant.  Additional  simulations  have  been  carried  out 
by  Amtec  Engineering,^  Yungster  at  the  University  of  Washington  and  later  at  NASA- 
Lewis^,  and  a  group  at  the  Naval  Research  Laboratory  (NRL).^®  Almost  all  researchers 
have  found  that  viscous  terms  must  be  included  in  the  simulations  and  that  turbulence 
modelling  affects  the  results.  While  the  fins  on  the  projectile  make  the  flowfield  three- 
dimensional,  most  numerical  studies  assume  axisymmetric  flow.  Despite  this  simplifica¬ 
tion,  good  agreement  with  experimental  data  has  been  obtained. 

Several  applications  of  the  concept  have  been  proposed.  One  of  the  best  uses  for  the 


3.  See  Knowlen,  Higgins,  and  Bruckner  (1994). 

4.  Recent  publications  are  Kruczynski  (1991)  and  Kruczynski,  Liberatore,  and  Kiwan  (1994). 

5.  See  Srulijes,  Smeets,  and  Seiler  (1992). 

6.  See  the  series  of  references  by  Nusca  (1991),  (1993)  and  (1994),  and  Kruczynski  and  Nusca 
(1992). 

7.  See  Sinha,  York,  Dash,  Drabczuk,  and  Rolader  (1992). 

8.  See  Soetrisno  and  Imlay  (1991)  and  Soetrisno,  Imlay  and  Roberts  (1992)  and  (1993). 

9.  See  Yungster,  Eberhardt,  and  Bruckner  (1991),  Yungster  (1991)  and  (1992),  and  Yungster  and 
Rabinowitz  (1993). 

10.  There  are  a  series  of  papers  by  Li,  Kailasanath,  and  Oran,  et.  al.,  in  the  list  of  references. 


3 


ram-accelerator  is  as  a  research  tool.  As  a  means  of  accelerating  masses  up  to  10  km/sec, 
the  device  can  provide  experimental  hypersonic  flow  data,  supplementing  ballistic  ranges, 
shock  tubes,  and  wind  tunnels.  Since,  by  its  nature,  it  relies  on  shock-induced,  and  at 
transdetonative  and  superdetonative  speeds,  supersonic  combustion,  the  ram-accelerator  is 
a  valuable  testbed  for  fundamental  research  of  those  phenomena.  Research  applications 
are  not  limited  to  gas  dynamics;  materials  analyses  and  impact  studies  are  also  plentiful.  A 
key  advantage  of  the  concept  is  that  scaling  to  large  projectile  diameters  and  masses  is 
straightforward.  Other  hypervelocity  acceleration  concepts  suffer  from  either  power  or 
material  property  limitations.  Because  scaling  is  not  an  issue,  the  ram-accelerator  is  being 

considered  as  an  economical  method  of  placing  acceleration-insensitive  payloads  into 
orbit. 

1.3  Unstarts 

The  establishment  of  a  steady,  thrust-producing  flowfield  upon  entrance  of  the  accel¬ 
erator  tube  IS  a  significant  achievement.^’  All  the  intricacies  of  the  starting  process  are  not 
yet  understood.  Typically  experimental  facilities  use  a  proven  ‘starting’  mixture  in  the  first 
tube  section  to  ensure  a  successful  injection  and  to  isolate  the  starting  transients  from  the 

remainder  of  the  tube.’^  The  dynamics  of  the  projectile  entering  the  accelerator  are  not 
considered  here. 

Once  a  stable  flowfield  is  attained,  steady  operation  may  be  interrupted  for  a  number 
of  reasons.  The  deterioration  of  the  cycle  and  consequent  loss  of  thrust  is  referred  to  as 
unstart.  Structural  failure  of  the  projectile  is  one  reason  for  unstarts.  Recent  efforts  at  the 
University  of  Washington’^  have  demonstrated  that  aluminum  forebodies  can  not 
withstand  the  high  heat  transfer  of  some  gas  compositions.  Operation  may  also  be 
inhibited  by  projectile  canting;  if  the  axes  of  the  accelerator  tube  and  projectile  are  not 


11.  Some  difficulties  experienced  at  ARL  are  described  in  Kruczynski,  Liberatore,  and  Kiwan 
(1994). 

12.  See  Knowlen,  Higgins,  and  Bruckner  (1994). 

13.  Experimental  documentation  of  structural  failures  is  given  in  Knowlen,  Higgins,  and  Bruckner 
(1994).  A  numerical  simulation  of  the  heating  of  a  projectile  forebody  is  presented  in  Chew  and 
Bruckner  (1994). 


aligned,  variations  in  the  shock  strength  from  the  projectile  tip  can  ignite  the  mixture 
prematurely.  Projectiles  with  five  fins  have  shown  better  performance  than  four-fin  types 
because  they  do  not  allow  as  much  canting.  Numerical  simulations  confirm  the 
sensitivity  of  stable  operation  to  slight  errors  in  projectile  alignment. 

Unstarts  can  result  even  for  well-aligned,  structurally  sound  projectiles.  Such  gasdy- 
namic  unstarts  generally  fall  into  one  of  three  categories,  all  of  which  end  with  significant 
energy  release  ahead  of  the  projectile.  Because  the  energy  is  released  prematurely  high 
pressures  are  experienced  on  the  forebody,  and  thrust  is  lost.  In  the  first  mechanism,  a  det¬ 
onation  wave  forms  in  the  wake  and  propagates  upstream  past  the  projectile  body.  It  either 
remains  attached  to  the  projectile  nose  or  propagates  through  the  tube  ahead  of  the  projec¬ 
tile.  Behind  this  detonation  wave  all  the  fuel  is  burned,  so  the  heat  release  occurs  too  early. 
This  type  of  unstart  generally  occurs  when  the  flow  in  the  wake  is  thermally  choked. 

In  the  second  mechanism,  the  projectile  boundary  layer  is  heated  by  the  impingement 
of  a  shockwave  on  the  afterbody.  Ignition  occurs  at  the  impingement  point  and  combus¬ 
tion  propagates  upstream  through  the  boundary  layer  to  the  forebody.  As  combustion 
moves  along  the  projectile,  it  spreads  to  the  rest  of  the  flowfield.  The  third  class  of  gasdy- 
namic  unstarts  sets  an  upper  speed  limit  on  the  projectile.  Above  some  speed  the  oblique 
shockwave  formed  at  the  nose  is  strong  enough  to  ignite  the  flow.  When  this  happens,  the 
energy  release  may  strengthen  the  shock  and  form  a  detonation  wave  that  propagates 
upstream.  Alternatively  the  energy  release  raises  pressure  on  the  forebody  enough  to 
decelerate  the  projectile.  It  is  important  to  note  that  combustion  in  the  forebody  boundary 
layer  or  behind  the  initial  shockwave  does  not  necessarily  result  in  an  unstart.  If  the  energy 
release  is  not  too  great,  the  rest  of  the  flowfield  will  not  ignite  until  further  downstream. 
Ideally,  however,  there  is  no  combustion  ahead  of  the  projectile. 

For  the  second  and  third  unstart  mechanisms  described,  projectile  geometry  modifica¬ 
tions  have  been  proposed  to  delay  or  eliminate  premature  combustion.  Rom  and  cowork- 
ers^^  propose  a  forward  facing  step  be  included  at  the  projectile  shoulder.  This  step  causes 

14.  See  Hinkey,  Burnham  and  Bruckner  (1992). 

15.  See  Soetrisno,  Imlay,  and  Roberts  (1993). 

16.  This  mechanism  can  be  seen  in  the  numerical  simulations  by  Yungster  (1991). 

17.  See  Rom  and  Avital  (1992),  and  Tivanov  and  Rom  (1993). 
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Figure  1.2  Double-cone  projectile  geometry. 

a  strong  shock  to  form,  which  ignites  the  gas  at  the  desired  location.  A  similar  variation 
offered  by  Yungster^^  consists  of  a  double-cone  forebody,  shown  in  Figure  1.2.  (The  fins 
are  not  shown  in  the  figure.)  The  first  cone  angle  is  small  to  weaken  the  initial  shock  and 
reduce  drag.  At  some  point  downstream  the  cone  angle  steepens,  forming  the  second  cone. 
The  shock  formed  by  the  second  cone  is  strong  enough  to  ignite  the  gas  mixture. 

1.4  The  Scope  of  Present  Work 

Ideally  the  numerical  simulation  of  ram-accelerator  flowfields  would  include  models 
for  all  the  physical  phenomena  that  occur.  A  detailed  chemical  model  including  a  multi¬ 
tude  of  species  with  reliable,  experimental  validation  at  ram-accelerator  conditions  would 
be  employed.  The  effects  of  intermolecular  forces  would  be  included  because  of  the  high 
pressures  encountered.  Turbulence  would  also  be  accurately  modelled.  The  physical 
domain,  including  the  wake,  would  be  represented  using  a  fine  mesh  for  the  entire  flow- 
field.  The  simulation  would  be  fully  three-dimensional  and  time-accurate. 

Unfortunately,  computer  speed  and  memory  are  limited;  some  models  are  too  expen¬ 
sive  computationally.  Others  are  of  limited  or  unverified  accuracy.  As  a  result  compro¬ 
mises  are  made  to  reduce  the  complexity  and  cost  of  the  simulations.  Despite  the 
compromises  reasonably  accurate  flowfield  solutions  have  been  obtained  by  many 
researchers.  Comparison  with  experimental  data  demonstrates  that  the  salient  physics  are 
included  in  the  models  and  important  trends  can  be  predicted. 

18.  See  Yungster  (1993). 
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In  this  thesis  ram-accelerator  fiowfields  are  simulated  over  a  range  of  velocities.  The 
objective  is  an  assessment  of  the  unstart  characteristics  of  the  double-cone  geometry.  In 
particular,  interest  is  centered  on  the  unstart  mechanisms  involving  combustion  in  the 
forebody  boundary  layer  or  ignition  by  the  initial  shockwave.  The  simulations  are  axisym- 
metric,  laminar,  steady,  and  do  not  encompass  the  projectile  wake.  Including  the  wake  dra¬ 
matically  increases  the  cost  of  the  simulations  but  does  not  play  a  significant  role  in  the 
unstart  mechanisms  to  be  studied.  For  similar  reasons  fins  on  the  afterbody  of  the  projec¬ 
tile  are  omitted.  The  gas  mixture  is  assumed  to  be  thermally  perfect,  reacting  methane-air. 

The  solution  advancement  scheme  has  been  used  for  reentry  flows,  but  not  for  com¬ 
busting  flows.  A  secondary  objective  is  to  further  establish  the  ‘modified  diagonal 
implicit’  scheme.^^  The  simulations  show  steady-state  solutions  for  two  reasons. 
Unsteady  simulations  demand  more  CPU  time  and  more  memory  for  storage  of  intermedi¬ 
ate  solutions.  The  diagonal  implicit  scheme  is  not  time  accurate  and  would  have  to  be 
replaced;  further  validation  would  not  be  possible.  The  details  of  the  physical  models,  the 
numerical  solution  technique,  and  the  resulting  fiowfields  are  discussed  in  the  following 
chapters. 


19.  The  modified  diagonal  implicit  scheme  is  presented  in  Candler  and  Olynick  (1991)  and  Hassan, 
Candler,  and  Olynick  (1992),  as  well  as  in  Chapter  HI. 
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GOVERNING  EQUATIONS 


11.1  Introduction 

In  this  chapter,  the  equations  governing  chemically  reacting  flowfields  are  introduced. 
The  assumptions  in  the  equations  are  also  discussed.  Some  assumptions  are  inherent  in  the 
equations  themselves,  while  in  other  cases  simplifications  are  made  to  make  solution  fea¬ 
sible. 

The  coupled,  nonlinear,  partial  differential  equations  that  describe  the  conservation  of 
mass,  momentum,  and  energy  are  collectively  called  the  Navier-Stokes  equations.  They 
are  augmented  for  application  to  chemically  reacting  flows.  The  definition  of  the  total 
energy  and  a  gas  law  specify  the  state  of  the  gas.  Constitutive  relations  provide  expres¬ 
sions  for  diffusion  of  mass,  momentum,  and  energy.  The  evaluation  of  empirical  coeffi¬ 
cients  required  for  these  expressions  is  explained  in  the  section  on  Transport  Properties. 
The  models  for  mass  transfer  between  chemical  species  are  described  in  the  Chemical 
Source  Terms  section.  Finally,  the  boundary  conditions  for  ram-accelerator  flowfields  are 
examined. 

11.2  Conservation  Equations  for  Reacting  Flow 

The  spatial  and  temporal  evolution  of  fluid  properties  is  described  by  the 
Navier-Stokes  equations.  This  set  of  nonlinear,  partial  differential  equations  specifies  the 
conservation  of  mass,  momentum,  and  energy.  The  fundamental  assumption  in  the  equa¬ 
tions  is  that  a  continuum  description  of  the  gas  is  applicable.  For  flows  in  which  the  char¬ 
acteristic  length  scale  is  similar  to  a  mean  free  path,  a  particle  description  is  appropriate. 
For  ram-accelerator  flowfields,  and  combustion  flowfields  in  general,  gas  densities  are 
high  enough  to  apply  the  Navier-Stokes  set.  By  the  same  reasoning,  no-slip  conditions  are 
valid  at  gas-solid  interfaces,  and  linear  constitutive  relations  are  accurate. 

The  Navier-Stokes  equations  can  be  augmented  to  model  chemically  reacting  flows. 
To  do  so,  the  conservation  equations  are  applied  to  each  chemical  species  in  a 
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multi-component  mixture  and  appropriate  source  terms  are  added  to  account  for  mass, 
momentum,  and  energy  exchanges  between  species.  The  Navier-Stokes  equation  set  can 
also  be  extended  to  model  thermal  nonequilibrium,  the  effects  of  an  electric  field,  and 
radiation,  but  these  effects  are  not  included  in  the  present  work. 

At  this  point  the  Navier-Stokes  equations,  expressed  in  index  notation,  are  introduced. 
The  natural  extension  of  index  notation  is  to  Cartesian  systems  although  it  can  be  applied 
to  any  orthogonal  coordinate  system.  In  Chapter  in  notation  specific  to  two-dimensional 
and  axisymmetric  geometries  will  be  addressed.  The  conservation  of  mass  for  species  ^  is 
expressed  as 


Tt 

The  conservation  of  momentum  for  species  s  is  given  by 

and  the  conservation  of  energy  is  given  by 


(2.2.1) 


(2.2.2) 


(2.2.3) 


In  theory  the  above  equations  could  be  solved  for  the  species  densities,  velocities,  and 
energies  to  obtain  flowfield  solutions.  Clearly  the  size  of  the  system,  and  thus  the  compu¬ 
tational  expense  of  the  solution,  grows  rapidly  with  the  number  of  species.  In  practiee 
some  simplifying  assumptions  are  made  to  reduce  the  number  of  equations. 

For  a  gas  model  comprised  of  ns  chemical  species,  there  are  ns  mass  conservation 
equations.  For  species  5,  the  density  is  denoted  by  and  the  velocity  in  the  j  direction 
is  u^  j .  The  net  production  or  destruction  due  to  chemical  reactions  is  ;  the  treatment  of 
is  elucidated  in  Section  II.6.  The  species  velocities  can  be  rewritten  as  the  mass-aver¬ 
aged  velocity,  Uj ,  plus  a  perturbation  called  the  diffusion  velocity,  j : 

\j  =  +  (2.2.4) 

The  diffusion  velocity  is  assumed  to  be  small  compared  to  the  mass-averaged  velocity, 
which  is  defined  by 
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(2.2.5) 


# 


—  ^c  i 

O 

ns  ^ 

where  the  total  density  is  the  sum  of  the  species  densities, 

p  =  Sp.- 

ns 


(2.2.6) 


Substituting  the  mean  and  diffusion  velocities  for  j  in  the  species  continuity  equations, 


^p5  a  ,  ,  a  , 


(2.2.7) 


By  summing  all  the  species  continuity  equations,  the  conservation  of  total  mass  is  recov¬ 
ered.  This  holds  because  total  mass  is  not  created  or  destroyed  by  chemical  reactions, 

=  0,  (2.2.8) 


ns 


and  because  there  is  no  net  diffusive  flux. 


Ip.v,,;  =  0.  (2.2.9) 

ns 

To  reduce  the  number  of  equations,  the  species  momentum  equations  are  summed  in 
each  coordinate  direction.  Dalton’s  law  of  partial  pressures  states  the  mixture  pressure  is 
the  sum  of  the  species  pressures: 


P  =  2?,.  (2.2.10) 

ns 

The  shear  stress  for  the  mixture  is  also  the  sum  of  the  species  quantities: 


i,j 


(2.2.11) 


With  the  assumption  that  all  particle  collisions  are  elastic,  the  sum  of  the  momentum 
source  terms  is  zero: 


=  0.  (2.2.12) 

ns 

The  elastic  collision  assumption  is  consistent  with  the  use  of  the  ideal  gas  equation  of 
state,  which  is  discussed  in  the  next  section. 

With  these  relations,  the  mass-averaged  momentum  equations  are 
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(2.2.13) 


|(P»,)  +^(p«,a^+p6y  +  t.  .)  =  0. 

The  symbol  5.  ^  is  the  Kronecker  delta.  By  the  mass-averaging,  a  significant  reduction  in 
the  number  of  equations  is  accomplished.  The  evaluation  of  the  momentum  exchange 
terms  .  is  no  longer  required,  which  is  also  a  notable  simplification.  However,  the  dif¬ 
fusion  velocities  must  now  be  approximated  because  there  are  fewer  equations  than 
unknowns. 

The  species  energy  equations  are  mass-averaged  in  the  same  manner  as  the  momen¬ 
tum  equations.  Since  the  total  energy  is  conserved, 

=  0-  (2.2.14) 

ns 

The  heat  flux  vector  is  the  sum  of  the  species  contributions: 

«<,)  =  (2-2-15) 

ns 

The  total  energy  equation  is  now  written  as 

+  +  +  =  0-  (2.2.16) 

ns  J 

A  detailed  definition  of  the  total  energy  per  unit  mass,  pE  ^  is  given  in  the  next  section. 
The  evaluation  of  the  heat  flux  vector,  the  shear  stress  tensor,  and  the  diffusion  velocities 
are  described  in  Section  II.4. 

For  reasons  which  will  be  thoroughly  discussed  in  Chapter  HI,  some  of  the  species 
mass  conservation  equations  are  replaced  by  elemental  mass  conservation  equations. 
Since  the  species  equations  implicitly  conserve  elemental  masses,  as  well  as  the  total 
mass,  the  physics  described  by  the  equation  set  are  unchanged;  the  equations  have  merely 
been  recast  in  different  variables.  There  are  still  a  total  of  ns  mass  conservation  equations. 

The  first  ne  species  mass  conservation  equations  are  replaced  by  elemental  mass  con¬ 
servation  equations,  but  the  remaining  ns  -  ne  species  equations  are  retained.  When  chem¬ 
ical  reactions  occur,  the  elemental  masses  are  not  affected;  consequently,  there  are  no 
source  terms  for  the  elemental  continuity  equations: 
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(2.2.17) 


The  density  of  element  e,  pg ,  is  related  to  all  the  species  densities  by  the  linear  combi¬ 
nation 


e=l,  ...,ne,  (2.2.18) 

S  =  1 

where  ne  is  the  number  of  elements.  The  total  density  is  now  defined  by 

ne 

P=Xp.-  (2.2.19) 

e  =  I 

The  diffusion  velocity  for  element  e,  v^  j ,  is  related  to  the  species  diffusion  velocities  by 

ns 

PA,;  =  E  «e,5pA;  e=\,...,ne.  (2.2.20) 

5=1 

Given  the  ne  elemental  conservation  equations  and  the  ns-ne  species  conservation  equa¬ 
tions,  the  replaced  species  densities  can  be  obtained  by 

ne  ns 

P^=EP^,^P^+  X  P^.fPf  s=\,...,ne.  (2.2.21) 

t=  I  t  =  ne+l 

The  coefficients  ^  and  ^  are  functions  of  the  elemental  and  species  molecular 
weights.  The  determination  of  these  coefficients  is  detailed  in  Appendix  A. 

With  the  explicit  inclusion  of  the  elemental  conservation  equations,  the  conserved 
variables  are  the  elemental  densities,  p^,  e  =  I,  ...,ne,  the  retained  species  densities, 
p^,  5  =  ne+l, ns,  the  momentum  in  each  spatial  dimension,  pu.,  and  the  total 
energy,  pE.  In  practice  the  densities  of  the  replaced  species,  p^,  ^  =  1,  ...,ne,  are  also 
calculated  because  many  quantities  are  more  difficult  to  evaluate  in  terms  of  the  elemental 
densities.  For  this  reason  the  equations  in  this  chapter  are  expressed  in  the  full  set  of  spe¬ 
cies  densities. 

In  summary,  a  set  of  nonlinear  partial  differential  equations  describing  fluid  flow  is 
presented.  The  set  consists  of  ns  equations  for  the  chemical  composition,  a  momentum 
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# 


# 


# 


equation  for  each  coordinate  direction,  and  a  total  energy  equation.  The  governing  equa¬ 
tions  are  expressed  on  a  per-unit-volume  basis. 

Conservation  of  elemental  mass: 


Conservation  of  species  mass: 


e  =  ...,ne 


(2.2.22) 


8  3 

Conservation  of  momentum: 

3  3 

gj(P“,)  +5^<P",“j+p5w  + V  =  0 


(2.2.23) 


(2.2.24) 


Conservation  of  total  energy: 


y  L  J 


(2.2.25) 


II.3  Equation  of  State 


In  this  section  the  relations  between  the  conserved  variables  and  the  non-conserved, 
intensive  variables  p  and  T  are  presented.  Once  the  conserved  variables  have  been  calcu¬ 
lated,  the  temperature  is  determined  iteratively  from  the  definition  of  the  total  energy.  An 
iterative  solution  is  required  because  of  the  treatment  of  the  vibrational  energy.  The  pres¬ 
sure  is  then  determined  from  the  ideal  gas  law. 

In  ram-accelerator  flowfields,  pressures  can  reach  hundreds  or  even  thousands  of 
atmospheres.  For  such  high  pressures,  intermolecular  forces  may  have  a  significant  effect 
on  the  state  of  the  gas.  However,  the  fluid  is  assumed  to  be  a  mixture  of  thermally  perfect 
gases  for  simplicity.  For  such  a  mixture,  the  ideal  gas  law  relates  the  pressure  to  the  tem¬ 
perature,  T,  and  density: 


P  =  'LPs  = 

ns 


(2.3.1) 


where  R  is  the  universal  gas  constant  and  the  molecular  weight  for  a  species  is  denoted 
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The  total  energy  is  comprised  of  translational-rotational,  vibrational,  chemical,  and 
kinetic  components: 


ns 


ns 


=  S  PA^^+  X  PA,i+  X  pX+  X  P^“^  -  (2.3.2) 

S=1  s=l  J=1  J=1 

The  energy  conservation  equation  requires  the  specific  enthalpy,  ,  for  each  species.  This 
is  given  by 


^v,  ^v,s 


I  1°  p. 

+  + 

\s<nv  p 


(2.3.3) 


In  both  these  equations,  ^  is  the  constant- volume,  translational-rotational  specific 
heat.  It  is  a  function  of  the  molecular  weight  and  internal  structure  of  the  species.  Contri¬ 
butions  from  the  translational  and  rotational  modes  are  summed: 

^v,s  ~  ^v,  trans,  i rot,  i  '  (2.3.4) 


The  contribution  to  ^  from  the  fully  excited  translational  mode  is  the  same  for  all  spe- 
cies,  .  The  rotational  mode  is  also  assumed  to  be  fully  excited: 


c 


V,  rot,  S 


0  for  monatomics 

—  for  diatomics  and  linear  polyatomics 

3  R 

z  —  for  3  -D  poly  atomics 


(2.3.5) 


There  are  three  rotational  modes  for  rotational  energy  for  molecules,  each  of  which  con- 
*17? 

tributes  to  the  specific  heat.  Because  rotation  about  the  bond-axis  is  insignificant  com- 
pared  to  the  other  two  axes  for  diatomic  and  linear  polyatomic  molecules,  only  two  modes 
contribute  to  the  rotational  specific  heat.  Monatomics  have  no  internal  structure  and  thus 
store  no  rotational  energy.  Note  that  the  vibrational  modes  also  contribute  to  the  heat 
capacities  of  molecules,  but  this  is  accounted  for  in  the  vibrational  energy  term  discussed 
below. 

The  vibrational  energy  modes  are  partially  excited.  The  total  vibrational  energy  is  the 
sum  of  the  species  vibrational  energies.  Monatomic  species  have  no  vibrational  energy 
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and  are  listed  last.  The  first  nv  species  are  polyatomic  and  consequently  store  energy  in 
vibrational  modes.  Each  of  these  species  has  ndeg(s)  significant  vibrational  modes  with 
characteristic  temperature  0^,  ^  ^ .  For  each  species  the  energy  in  all  modes  is  summed 
with  degenerate  modes  counted  multiple  times: 


E..  .  = 


ndeg  (^) 

—  y 

^  m=l 


V,  s,  m 


exp 


0 


V,  Sy  m 


(2.3.6) 


Note  that  while  the  vibrational  modes  are  partially  excited,  they  are  in  equilibrium  with 
the  translational  and  rotational  energies.  As  a  result  a  single  temperature  characterizes  the 
internal  energy  of  the  gas  particles,  and  a  separate  partial  differential  equation  for  the 
vibrational  energy  is  not  required.  The  temperature  must  be  found  iteratively  from  the  def¬ 
inition  of  the  energy;  Equation  (2.3.2)  can  not  be  manipulated  to  a  closed  form  for  T 
because  of  the  vibrational  energy  terms.  Finally,  curve  fits  for  thermodynamic  properties 
are  not  required  with  this  formulation  of  the  total  energy;  the  0^  ‘s  specify  the  variance  of 
the  heat  capacities  with  temperature. 

The  third  term  of  Equation  (2.3.2)  establishes  a  reference  level  of  energy.  When  a  spe¬ 
cies  is  produced  or  consumed  by  chemical  reaction,  the  energy  released  or  absorbed  is 
accounted  for  by  the  enthalpy  of  formation,  .  Kinetic  energy  due  to  the  macroscopic 
motion  of  the  gas  is  included  in  the  fourth  term.  The  enthalpies  of  formation  are  found  in 
the  literature  and  are  listed  for  each  species  in  Appendix  B,  as  are  the  molecular  weights 
and  characteristic  vibrational  temperatures.^ 


II.4  Constitutive  Relations 

Constitutive  relations  give  expressions  for  the  diffusion  velocities,  the  stress  tensor, 
and  the  heat  flux  vector.  The  diffusion  of  mass,  momentum,  and  energy  are  dependent  on 
first  order  gradients  of,  respectively,  concentration,  velocity,  and  temperature. 

In  this  work  mass  diffusion  results  from  flowfield  concentration  gradients.  Pressure 


1.  Data  is  from  Stull  (1965)  and  (1971),  Dixon-Lewis  (1984),  and  McBride,  Heimel,  Ehlers,  and 
Gordon  (1963). 


and  temperature  driven  mass  diffusion  are  neglected.  Additionally  it  is  assumed  that  diffu¬ 
sion  is  binary;  that  is,  species  ^  diffuses  into  a  bath  of  similar  particles.  The  binary  diffu¬ 
sion  assumption  is  reasonable  when  the  molecular  weights  of  all  species  are  about  the 
same.  With  these  assumptions,  the  diffusion  velocities  are  found  from  Pick’s  Law 


(2.4.1) 


where  is  the  binary  diffusion  coefficient  and  is  the  mass  fraction  of  species 


X  = 


(2.4.2) 


More  accurate  methods  for  the  determination  of  the  diffusion  velocities  are  available  but 
are  computationally  more  intensive.^ 

Shear  stresses  in  the  fluid  are  a  consequence  of  velocity  gradients.  The  gas  mixture  is 
treated  as  a  Newtonian  fluid,  for  which 


du-  du.  I  du. 

The  Stokes  hypothesis  for  the  bulk  viscosity  is  used  to  specify  X : 


(2.4.3) 


(2.4.4) 


The  heat  flux  vector  arises  from  temperature  gradients  in  the  flow.  Fourier’s  Law  stipu- 


ar 


(2.4.5) 


Evaluation  of  the  binary  diffusion  coefficient,  ,  the  coefficient  of  molecular  viscosity, 
(1 ,  and  the  heat  transfer  coefficient,  K ,  is  discussed  in  the  next  section. 

11.5  Transport  Properties 

The  coefficients  D^,  p,  and  K  are  required  for  the  calculation  of  the  viscous 


2.  A  method  involving  a  series  solution  for  the  diffusion  velocities  is  outlined  in  Oran  and  Boris 
(1987),  pp.  253-255,  and  is  also  described  in  Wilson  (1991). 
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phenomena  described  by  the  constitutive  relations.  Kinetic  theory  provides  some 
information,  but  ultimately  experimental  measurements  are  necessary  to  arrive  at  values 
for  the  coefficients. 

Attention  is  first  given  to  the  coefficient  of  viscosity.  The  approach  taken  consists  of 
two  parts.  A  first-order  approximation  to  the  coefficient  is  obtained  for  each  species  as  a 
pure  gas.  Then  these  species  quantities  are  combined  using  a  mixing  rule  to  get  a  value  of 
p,  for  the  mixture. 

From  kinetic  theory. 


2.6693x10  ^ 


(2,2)* 

ss 


(PN) 


(2.5.1) 


The  subscript  1  specifies  the  species  viscosity  coefficient  is  a  first-order  approximation.  In 
this  equation,  T*  is  the  reduced  temperature,  defined  by 

kj,T 

=  —  (2.5.2) 

e. 


where  kg  is  the  Boltzmann  constant.  The  constants  and  are  Lennard-Jones  potential 

(2  2)* 

parameters.  The  symbol  ’  is  the  elastic  collision  integral; 

^(2,2)*  ^  ^^^(1,1)*^  (2.5.3) 
(1  1)* 

For  each  species.  A*  and  Q  ’  are  calculated  from  polynomials  in  T* .  The  Len- 
nard-Jones  parameters  and  curve  fit  coefficients  are  taken  from  Dixon-Lewis^  and  are 
repeated  in  Appendix  B.  With  the  listed  values,  the  units  of  [u.  ] ,  are  kg/m-sec. 

Wilke‘s  semi-empirical  mixing  rule^  is  employed  to  obtain  a  mixture  viscosity  coeffi¬ 
cient  from  the  species  coefficients.  The  rule  is 


11  = 


ns 


1 


ns 


Lr=  1 


(2.5.4) 


3.  See  Dixon-Lewis  (1984). 

4.  See  Wilke  (1950). 
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where 


mJ  ((^^]J  [mJ 


(2.5.5) 


The  mole  fraction  of  species  s  is 


^  p,V' 


7  =  —  V  — 
^  M,  ^  M 


(2.5.6) 


The  mixture  heat  transfer  coefficient  is  also  found  by  first  determining  pure  species 
values  and  then  combining  them  with  the  Wilke  rule.  The  species  heat  transfer  is  related  to 
the  viscosity.  It  is  assumed  that  the  transport  of  translational  energy  can  be  correlated  with 
the  velocity,  but  for  internal  energy  there  is  no  such  correlation.  This  assumption  is  called 


the  Eucken  correlation.^  The  result  is 


’-c  +  r  4-  c 

2  V,  trans,  s  v,  rot,  s  v,  vib,  s  j  * 


(2.5.7) 


Previously  ^  and  ^  have  been  discussed.  The  vibrational  specific  heat  is  just 


j  .  .  V,  5,  /n  V,  5,  m 


'v,vlb,s  Jj.  hi 


V,  5,  m  1 


(2.5.8) 


As  with  the  viscosity,  Wilke’s  rule  is  used  to  find  the  mixture  heat  transfer  coefficient: 


K  =  y  — — 

Zmj  ns 

Spa,, 


(2.5.9) 


The  binary  diffusion  coefficient  is  estimated  by  assuming  a  constant  Schmidt  number 


of  1/2: 


=  -^  =  I. 

pD  2 


(2.5.  lO) 


5.  See  Vincenti  and  Kruger  (1965),  pp.  19-21. 
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where  it  is  further  assumed  that  all  species  diffuse  at  the  same  rate  and  have  the  same  D. 
Because  molecular  and  diatomic  hydrogen  have  molecular  weights  much  lower  than  the 
other  chemical  species  in  the  gas  model,  they  should  diffuse  more  rapidly.  Some  investi¬ 
gators  use  higher  values  for  and  to  account  for  the  inaccuracy  of  the  binary  diffu¬ 
sion  assumption,  which  is  the  source  of  the  underestimated  diffusion  for  the  less  massive 
species.  However,  this  practice  sacrifices  the  conservation  of  mass  and  is  not  implemented 
here. 

II.6  Chemical  Source  Terms 

The  chemical  source  terms  account  for  the  change  in  mixture  composition  due  to 
chemical  reaction.  A  nonequilibrium  treatment  is  appropriate  because  the  characteristic 
chemical  time  scale  is  on  the  same  order  as  the  characteristic  fluid  time  scale.  An  assumed 
set  of  chemical  reactions  is  used  to  model  the  change  in  species  mass  fractions  throughout 
the  flow.  There  are  several  drawbacks  to  the  model,  but  these  are  outweighed  by  the  sav¬ 
ings  in  computational  resources. 

A  general  chemical  reaction  can  described  by  the  stochiometric  equation 

n  n 

(2.6.1) 

i  = i  i  =  I 

In  this  equation,  species  i  has  the  chemical  symbol  C. .  The  stochiometric  coefficients  for 
the  reactants  are  denoted  by  v/  and  for  the  products  by  v /' .  There  are  a  total  of  n  chemi¬ 
cal  species  involved  as  reactants,  products,  or  collision  partners.  The  general  equation 
above  in  fact  represents  two  opposing  reactions.  The  forward  reaction  proceeds  from  reac¬ 
tants,  on  the  left,  to  products,  on  the  right.  For  the  backward  reaction  the  reactants  and 
products  are  reversed,  and  the  reaction  proceeds  from  right  to  left.  The  notation  can  be 
used  for  a  set  of  opposing  reactions.  In  that  case, 

ns  ns 

^  ’  (2.6.2) 

i  =  1  i=l 

where  r  refers  to  the  particular  reaction.  When  a  set  of  reactants  is  considered,  n  is  the 
number  of  chemical  species  in  the  gas  model  so  n  =  ns.  If  species  J  does  not  participate  in 
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reaction then V .  '  =  v.  "  =  0. 

’  J,S  J,S 

According  to  the  law  of  mass  action,  the  reaction  rate,  RR^,  is  given  by 

ns  ,  ns 

RR.  =  h.rU.  [C,l  [C'/'"  -  (2-6-3) 

*  =  1  z  =  1 

where  ^  and  ^  are  respectively  the  forward  and  backward  reaction  rate  constants  for 
reaction  r.  The  square  brackets  denote  the  concentration  (particles/volume)  of  the  enclosed 
species.  The  reaction  rate  and  the  chemical  source  terms  are  closely  related;  for  species  s 
the  net  production  or  consumption  of  mass,  after  all  nr  chemical  reactions  are  considered, 
is 


^  r=  1 


(2.6.4) 


The  forward  rate  constant  is  expressed  in  the  modified  Arrhenius  form. 


k^^^  =  A^r^'exp[  - 


0 


d,  r 


(2.6.5) 


The  characteristic  temperature  of  dissociation,  0^  the  temperature  exponent,  rj^,  and 
the  pre-exponential  factor,  A^,  are  parameters  determined  by  experimental  methods  and 
are  found  in  the  literature. 

The  backward  reaction  rate  constant  is  determined  from  the  forward  rate  constant  and 
the  equilibrium  constant  based  on  concentration,  ^ : 


K 


eq,  r 


-  hi 

b,r 


(2.6.6) 


The  equilibrium  constant  is  calculated  using  the  Gibbs’  free  energies  of  each  species  that 
participates  in  the  reaction.  From  chemical  thermodynamics, 


^eq,.=  (/?70'^^exp 


f 


AG 

I 

'rt 


(2.6.7) 


where 


z  =  1  z  =  1 


(2.6.8) 
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and 


ns 


i  =  1 


i  =  1 


(2.6.9) 


The  molar  standard  Gibbs’  free  energy  of  formation  for  a  species,  g°  ,  is  calculated  from  a 
polynomial  in  T.  The  curve  fit  coefficients  are  determined  experimentally^  and  are  listed  in 
Appendix  B. 

The  assumed  set  of  chemical  reactions  is  called  a  chemical  kinetics  mechanism.  In 
practice,  a  particular  mechanism  is  valid  for  a  given  range  of  conditions  or  flow  situations. 
Species  and  reactions  are  included  or  omitted  based  on  their  significance  in  the  mecha¬ 
nism’s  intended  application.  The  significance  of  a  particular  reaction  or  species  is  not 
always  clear. 

In  addition  to  the  reactions,  a  mechanism  specifies  the  parameters  for  determining  the 
forward  rate  constants.  The  parameters  are  chosen  to  match  relevant  experimental  data. 
Although  the  modified  Arrhenius  form  for  the  rate  constants  can  be  derived  from  kinetic 
theory,  the  actual  determination  of  A^,  T]^,  and  6^  ^  is  difficult.  Because  of  uncertainties 
in  these  rate  parameters  and  the  closely  related  issue  of  the  importance  of  particular  reac¬ 
tions,  many  different  mechanisms  have  been  proposed  for  hydrocarbon  combustion. 

For  the  present  studies,  methane  is  used  as  the  fuel  and  oxygen  as  the  oxidizer.  Nitro¬ 
gen  is  also  included  in  the  mixture  as  a  diluent.  Much  work  has  been  done  to  determine  the 
behavior  of  hydrocarbon  fuels  in  internal  combustion  and  gas  turbine  engines.  Unfortu¬ 
nately  the  resulting  kinetics  mechanisms  are  complex,  including  as  many  as  50-i-  species 
and  hundreds  of  reactions.  Since  the  required  computational  time  and  memory  scale 
roughly  as  the  square  of  the  number  of  species,  complex  mechanisms  can  not  be  used 
except  for  very  simple  flows.  Furthermore,  most  combustion  mechanisms  have  been 
developed  for  modeling  subsonie  flames  at  or  near  atmospheric  pressures.  Ram-accelera¬ 
tor  flowfields  are  supersonic  or  hypersonic,  and  pressures  range  to  hundreds  or  even  thou¬ 
sands  of  atmospheres. 

A  quasi-global  model  is  used  in  the  current  research.  In  the  quasi-global  methane 


6.  The  Gibbs’  free  energy  data  is  listed  in  McBride,  Heimel,  Ehlers,  and  Gordon  (1963). 
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reaction  mechanism  of  Westbrook  and  DryerJ  methane  and  oxygen  form  carbon 
monoxide  and  hydrogen  in  a  global  reaction: 

CH^  +  ^0^<^C0  +  2H^. 

A  set  of  elementary  reactions  models  the  oxidation  of  CO  and  H2  ■  The  global  reaction 
represents  a  number  of  elementary  reactions  in  a  single  step.  It  is  fundamentally  different 
from  the  elementary  reactions  it  approximates  because,  in  reality,  CH^  and  Oj  do  not 
form  carbon  monoxide  and  hydrogen  directly.  In  an  elementary  reaction,  the  reactants  are 
converted  to  products  without  intermediate  steps. 

The  quasi-global  mechanism  involves  the  reacting  species  CH^,  CO,  CO 2,  O2, 
H2O,  H2,  OH,  O,  and  H.  Nitrogen,  N2,  is  also  included  as  a  diluent  in  the  gas  model. 
The  mechanism  includes  12  reactions,  including  the  global  reaction.  The  reactions  are 

CO  +  O  +  M  CO 2  "I"  M 

0H  +  M<^0  +  H  +  M 
O2  +  M0O  +  O  +  A/ 

H2  +  M<^H  +  H  +  M 
H20  +  M<^0H  +  H  +  M 
CH^  +  ^^2  <=>  CO  +  2H2 
C0  +  0H^C02  +  H 
CO  +  O2  ^ 

H+02<^0H+0 
0  +  H2^0H  +  H 
0H  +  H2^H  +  H20 
H20  +  0^0H+0H 

M  is  a  collision  partner;  it  can  be  any  species.  The  rate  parameters  are  listed  in  Appendix 
B.  The  sixth  reaction  is  the  global  reaction.  The  original  reference  also  includes  the  radi¬ 
cals  HO 2  and  H2O2 .  These  species  are  important  for  modelling  ignition  at  lower  temper¬ 
atures.  Ignition  in  ram-accelerators  occurs  at  higher  temperatures;  additionally,  the  use  of 
the  quasi-global  model  nullifies  their  effects  on  ignition  in  this  mechanism,  as  explained 


7.  See  Westbrook  and  Dryer  (1981). 


below.  Therefore,  HO2  and  H2O2  and  the  reactions  they  participate  in  are  neglected  in 

o 

this  thesis. 

The  global  reaction  has  a  phenomenological  forward  reaction  rate,  given  by 

=  A6exp(-  ^-^^[CH^]\02]\  (2.6.10) 

The  form  of  the  rate  is  the  same  as  the  modified  Arrhenius  expression,  except  the  concen¬ 
tration  exponents  a  and  b  are  not  the  stochiometric  coefficients;  they  are  parameters  cho¬ 
sen  to  match  experimental  data.  In  this  case  a  =  -0.3  and  b  =  1.3.  Since  the  fuel 
exponent  is  negative,  the  forward  reaction  rate  approaches  infinity  as  the  fuel  concentra¬ 
tion  approaches  zero.  This  can  cause  problems  for  numerical  schemes.^  To  alleviate  the 
difficulty  Westbrook  and  Dryer  suggested^®  the  opposing  reaction  be  included  with  the 
rate 

RR6,»  =  .  (2.6.11) 

In  the  backward  reaction,  a'  =  -0.5  and  b'  =  1.3 .  For  the  global  reaction  the  reaction 
rate  is  thus 

rr6  =  RR6./-RR6.fc.  (2.6.12) 

and  the  methane  concentration  is  bounded. 

Because  several  elementary  reactions  are  approximated  with  a  single  global  reaction, 
the  quasi-global  mechanism  is  not  as  accurate  as  a  detailed  mechanism.  In  the  Westbrook 
and  Dryer  model,  the  global  reaction  replaces  elementary  reactions  associated  with  chain 
initiation  and  hydrocarbon  radicals,  and  as  a  consequence  the  details  of  the  reaction  front 
are  lost.  In  particular,  the  temperature  and  species  mass  fractions  through  the  reaction 
front  will  not  be  correctly  predicted.  On  the  other  hand,  the  burned-gas  composition  and 
temperature  should  be  accurately  resolved  because  the  oxidation  of  carbon  monoxide  and 
hydrogen  are  modelled  by  elementary  reactions.  Additionally,  the  range  of  applicability  of 

8.  Soetrisno  and  Imlay  (1992)  used  the  quasi-global  model  for  ram-accelerator  studies.  They  omit¬ 
ted  the  HO 2  and  H2O2  species  and  the  associated  reactions  without  adverse  effects. 

9 .  See  Westbrook  and  Dryer  (1981). 

10.  The  opposing  reaction  and  rate  listed  here  are  proposed  in  Coffee  (1985). 


the  mechanism  is  further  reduced  by  representing  several  elementary  reactions  with  a  sin¬ 
gle  global  reaction.  The  authors  point  out  the  pre-exponential  factor  should  be  taken  as 
an  initial  estimate.  Established  data  at  conditions  near  those  to  be  simulated  should  be 
used  to  calibrate  .  The  details  of  such  a  calibration  are  discussed  in  Chapter  IV. 

It  has  also  been  reported  that  the  quasi-global  is  sensitive  to  the  modeling  of  the  vis- 
cous  terms.  The  flame  speed  depends  to  some  degree  on  the  diffusion  of  mass,  momen¬ 
tum,  and  energy,  so  the  initial  calibration  of  the  global  reaction  rate  relies  on  the  treatment 
of  the  viscous  terms  in  the  original  reference.  Since  the  fuel  and  oxidizer  are  premixed  in 
the  ram-accelerator  and  since  the  pre-exponential  coefficient  will  be  recalibrated  anyway, 
the  sensitivity  to  the  viscous  modeling  is  less  critical  in  the  present  work. 

The  advantage  of  the  quasi-global  mechanism  is  that  it  includes  only  9  reacting  spe¬ 
cies  and  12  reactions.  Computationally  it  is  much  less  expensive  than  a  full,  detailed 
mechanism.  The  disadvantages  of  the  quasi-global  model  have  been  noted;  it  is  an  approx¬ 
imation  to  the  behavior  a  detailed  mechanism  predicts  at  ‘design’  conditions.  However,  at 
ram-accelerator  conditions,  the  available  detailed  mechanisms  are  no  less  suspect  but 
much  more  expensive  than  the  quasi-global  model.  For  this  reason  the  quasi-global  mech¬ 
anism  is  employed  in  this  thesis. 

II.7  Boundary  Conditions 

In  the  laboratory,  the  ram-accelerator  projectile  travels  through  the  pressurized  accel¬ 
erator  tube.  The  tube  has  a  constant-diameter  and  cross-section  and  is  filled  with  a  homo¬ 
geneous  mixture  of  fuel,  oxidizer,  and  diluent.  Before  the  arrival  of  the  projectile  the  gas  is 
quiescent.  Different  sections  of  the  tube  are  divided  by  diaphragms  and  can  be  indepen¬ 
dently  pressurized  and  fueled.  To  obtain  the  desired  performance,  the  pressure  and  mix¬ 
ture  composition  are  varied  in  the  different  tube  sections.  The  acceleration  of  the  projectile 
and  the  flowfield  it  induces  are  clearly  transient  in  the  laboratory  frame  of  reference. 

The  simulation  considers  the  steady-state  flowfield  only;  it  is  assumed  that  the 
transient,  real-world  flowfield  is  accurately  approximated  by  a  series  of  steady-state. 


11.  See  Westbrook  and  Dryer  (1981). 

12.  See  Coffee  (1985). 
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numerical  solutions.  In  the  simulation,  the  projectile  is  still  and  the  gas  flows  past  it. 
Ahead  of  the  projectile,  the  gas  velocity  is  constant  and  equal  to  the  projectile  speed  being 
simulated.  At  the  projectile  surface,  the  continuum  assumption  stipulates  the  gas  velocity 
is  zero  with  respect  to  the  projectile.  Similarly  the  no-slip  condition  applies  at  the 
accelerator  tube  wall.  In  the  simulation,  the  wall  moves  at  the  same  speed  as  the 
free-stream  gas  because  in  the  laboratory  the  projectile  moves  from  one  end  of  the  tube  to 
the  other. 

As  the  projectile  accelerates,  it  is  heated  by  skin  friction  and  the  high  temperatures  of 
the  burned  gas.  The  dynamics  of  the  temperature  field  within  the  projectile  are  not  consid¬ 
ered  in  this  study,  and  until  reasonable  estimates  of  the  wall  and  projectile  surface  temper¬ 
atures  are  available,  adiabatic  flowfield  boundaries  are  assumed.  For  simplicity  a 
non-catalytic  assumption  is  employed  for  the  tube  and  projectile  surfaces. 

With  the  above  boundary  conditions,  which  specify  composition,  velocity,  and  tem¬ 
perature  at  the  projectile  and  accelerator  tube  surfaces,  the  governing  equations  can  be 
solved  over  the  interior  of  the  computational  domain. 
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13.  Nusca  (1993)  points  out  that  the  wall  temperature  specification  significantly  affects  the  flow- 
field,  but  that  until  experimental  measurements  are  available  an  adiabatic  wall  condition  is 
appropriate. 
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CHAPTER  III 


NUMERICAL  METHOD 


111.1  Introduction 

The  equations  introduced  in  Chapter  n  are  solved  numerically.  The  partial  differential 
equations  are  rewritten  as  algebraic  difference  equations,  for  which  many  solution  tech¬ 
niques  have  been  established.  Spatial  terms  are  upwinded  or  centrally  differenced.  An 
implicit,  iterative  method  is  used  to  advance  the  flowfield  from  an  initial  condition  to  the 
steady-state  solution.  Evolution  of  the  flow  properties  is  driven  by  the  conditions  at  the 
boundaries  of  the  computational  domain. 

111.2  Vector  Form  of  the  Governing  Equations 

The  governing  equations  are  written  in  vector  form.  Some  changes  are  made  in  the 
notation.  Specifically,  the  indices  specifying  coordinate  directions  are  no  longer  used. 
With  the  modifications,  the  equations  describe  two-dimensional  or  axisymmetric  flow. 
The  partial  differential  equation  set  is  expressed  as 

=  W.  (3.2.1) 

dt  dx  y’'dy^  ^ 

For  two-dimensional  Cartesian  coordinates,  r  =  0. 

If  r  =  1 ,  the  equations  describe  axisymmetric  flow  in  cylindrical  coordinates.  Then 
the  geometry,  and  consequently  the  flow,  is  symmetric  about  some  centerline.  The  y-coor- 
dinate  is  the  radial  distance  from  this  centerline.  The  source  vector,  W,  contains  some  new 
terms  that  arise  from  the  coordinate  system.  Other  modifications  due  to  the  choice  of  the 
coordinates  will  be  discussed  when  necessary. 

The  conserved  variables  comprise  the  vector  U: 

U  =  (Pi, . . . ,  p„,,  . . . ,  p„,,  PM,  pv,  pF)  ^  .  (3.2.2) 

The  mass-averaged  velocities  in  the  x-  and  y-directions  are  u  and  v.  The  conserved  vari¬ 
ables  are  the  same  in  both  coordinate  systems. 


The  flux  vectors,  F  and  G,  are  separated  into  inviscid  and  viscous  components: 
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The  elemental  diffusion  velocities  in  the  jc-direction  have  symbol  ,  while  the  species 
diffusion  velocities  are  given  by  .  Similarly,  for  the  y-direction,  elemental  and  species 
diffusion  velocities  are  denoted  v  and  v  ,  and  the  flux  vector  is 
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(3.2.4) 


The  source  vector  contains  terms  that  account  for  mass  transfer  between  species  due  to 
chemical  reactions  and  a  term  that  arises  from  the  axisymmetric  coordinate  system: 

W  =  (0, . . . ,  0,  . . . ,  0,r(p-  Xqq)  /y,  0)  ^  .  (3.2.5) 


The  first  ne  entries  are  zero  since  the  elemental  mass  densities  are  unchanged  by  reactions. 
The  source  term  in  the  y-momentum  equation  results  from  the  axisymmetric  coordinates. 


III.3  Finite  Volume  Discretization 


The  partial  differential  vector  equation  (3.2.1)  can  be  discretized,  resulting  in  a  large 
system  of  algebraic  difference  equations.  A  finite-volume  formulation  over  the  computa¬ 
tional  mesh  is  used.  Inviscid,  viscous,  and  implicit  terms  arise  in  the  discretization  and  are 
handled  differently.  The  discussion  opens  with  the  definition  of  the  computational  domain. 

The  strategy  taken  by  the  numerical  method  is  to  represent  the  physical  domain  by  a 
limited  number  of  points.  At  those  points,  which  make  up  the  computational  domain,  the 
governing  equations  are  solved.  The  flowfield  properties  anywhere  in  the  physical  domain 
can  then  be  determined  by  interpolation  from  the  points  in  the  computational  domain.  The 
underlying  assumption  in  this  approach  is  that  the  physical  domain  can  be  accurately  rep¬ 
resented  by  this  collection  of  discrete  points.  If  enough  points  are  included,  the  assump¬ 
tion  is  valid. 

To  represent  complicated  geometries,  general  curvilinear  grids  are  preferred  over  rect- 
angular-celled  meshes.  The  equations  are  transformed  from  the  physical  space  (jc,  y)  to  the 
computational  space  (^,  T) ).  The  calculation  concerns  the  values  of  the  properties  at  the 
discrete  points  that  have  integer  values  of  ^  and  T) .  These  integer  values  are  given  by  the 
indices  i  and  j,  respectively.  A  generic  grid  cell  is  shown  in  Figure  3.1.  In  the  finite- vol¬ 
ume  formulation,  the  coordinates  x-  j  and  y.  j  refer  to  the  physical  location  of  the  lower 
left  corner  of  the  computational  cell.  However,  for  all  the  properties,  such  as  p .  . ,  , 

u.  J ,  and  so  forth,  the  indices  i  and  j  refer  to  values  at  the  cell  center. 

For  each  cell.  Equation  (3.2. 1)  is  integrated  to  give 

Hf-®  ^  J  7  ^  =  J  wn ,  (3.3.1) 

£2  0^  £2 

where  Q  is  the  volume  of  the  cell.  On  a  general  curvilinear  grid,  the  cells  are  not  aligned 
with  the  coordinate  axes.  For  this  reason  the  radius  included  for  the  axisymmetric  coordi¬ 
nates  is  applicable  to  both  the  F  and  G  fluxes. 

The  divergence  theorem  is  invoked  for  the  flux  terms: 

J  =  “r  J  [/Fi  +  yGj]  •  dS.  (3.3.2) 

£2  ^  ^5 
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Figure  3.1  Property  value  and  grid  definition  locations  for  a  representative  cell. 

The  volume  integral  is  recast  as  a  surface  integral,  where  the  surface,  5 ,  encloses  the  vol¬ 
ume  Q ,  and  i  and  j  are,  respectively,  unit  normal  vectors  in  the  jc-  and  y-directions.  In  the 
current  research  the  grid  is  fixed  throughout  the  computation,  so  Equation  (3.3.1)  can  be 
expressed  as 

J  +  yGj]  dS  =  aw.  (3.3.3) 

y  s 

The  temporal  term  is  forward-differenced,  and  the  surface  integral  is  replaced  as  a  sum; 
then 

^  ^  ^  X  j  ■  ^  (3-3.4) 

J  Sides 

where  Slf  =  -  if .  The  superscript  here  denotes  the  iteration  of  the  numerical 

scheme.  The  time  level  of  the  fluxes  and  the  source  vector  will  be  specified  below. 

At  this  point  consider  the  summation.  For  each  cell  side  the  appropriate  S  is  used.  For 
example,  consider  the  i  +  1/2  cell  side: 
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(3.3.5) 


‘^1+1/2  -  (  {  ^ix  ^  +  ^iy  h  -S'] 

^  h+\/2 

where  and  j  are  direction  cosines  and  S  is  the  area  of  the  cell  side.  Similarly  for  the 
j  +  1/2  cell  side, 


^j^i/2  =  {{^jx'^i  +  ^jy'j}s]  .  (3.3.6) 

The  direction  cosines  and  cell  side  areas  are  closely  related  to  the  mesh  metrics  and  their 
inverses.  The  calculation  of  these  quantities  is  discussed  in  Appendix  C,  along  with  other 
grid-related  variables. 

With  these  definitions,  fluxes  can  be  defined  in  curvilinear  coordinates.  Denoting  the 
rotated  fluxes  by  F'  and  G' , 

F'  =  Fs.;  +  Gs,;  (3.3.7) 

and 


G'  =  Fsj^  +  Gsjy'.  (3.3.8) 

The  definition  of  the  contravariant  velocities  depends  on  the  cell  face  being  considered. 
For  i  +  1/2  cell  faces, 

u'  =  us.^'  +  vs.y'  v'  =  -us.y'  +  vs.^',  (3.3.9) 

while  for  j  +  1/2  cell  faces, 

u'  =  usjy'-vsjj  v'  =  usj^  +  vsjy'.  (3.3.10) 

Defining  contravariant  velocities  simplifies  the  notation. 

Now  the  equation  to  be  solved  at  each  point  i,j  is 
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(3.3.11) 


where  the  time  level  n+l  has  been  specified  for  F' ,  G' ,  and  W.  The  cell  centroid  radius 
is  denoted  by  y^ ;  it  is  not  evaluated  at  the  same  location  as  the  radii  multiplying  the  fluxes 
F  and  G.  For  the  quantities  with  a  single  subscript,  the  other  index  is  implied;  it  is  either  i 
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or  j,  whichever  is  absent. 

Since  the  conserved  variables  at  iteration  n  +  1  are  not  yet  known,  the  fluxes  at  n  +  1 
are  estimated  based  on  the  solution  at  iteration  n.  The  fluxes  are  linearized  by 


F,n+\  ^  5f/*  = 

oU 


(3.3.12) 


G,n+i  ^  ^ 

aU 


(3.3.13) 


The  matrices  A  and  B  are  the  flux  Jacobians.  The  source  vector  can  be  linearized  in  the 


same  manner: 


+  ^  6lf  =  w”  +  z”6c/‘  , 


(3.3.14) 


where  Z  is  the  chemical  Jacobian. 

With  the  linearizations  above,  Equation  (3.3.11)  is  written 


yc,i,pi,j{ 


(3.3.15) 
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The  left-hand  side  of  this  equation  is  referred  to  as  the  implicit  side.  The  right-hand  side  is 
called  the  explicit  side,  and  is  defined  by 
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(3.3.16) 


The  evaluation  of  the  right-hand  side  is  discussed  in  the  next  section,  while  the  treatment 
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of  the  implicit  side  is  covered  in  Section  111.5.  A  technique  for  solving  Equation  (3.3.15)  is 
explained  in  Section  III.6. 

III.4  Evaluation  of  the  Explicit  Terms 

First,  consider  the  inviscid  fluxes  E'”  .^j/2  G'"  j  +  i/2-  For  supersonic  and 

hypersonic  flows  the  inviscid  terms  are  hyperbolic,  characterized  by  discontinuities  such 
as  shock  waves.  Upwind,  shock-capturing  methods  have  been  developed  to  accurately 
resolve  such  discontinuities.  The  original  upwind  methods  eliminated  spurious,  numerical 
oscillations  associated  with  central-difference  schemes  but  were  only  first-order  accurate. 
These  upwind  methods  have  been  extended  to  global,  second-order  accuracy  by  the  addi¬ 
tion  of  nonlinear  numerical  dissipation;  a  subset  of  these  modern  methods  contains  the 
Total  Variation  Diminishing,  or  TVD,  schemes.  For  the  simulations  presented  in  the  fol¬ 
lowing  chapters,  the  upwind  method  of  Roe^  is  extended  to  second-order  accuracy  by  the 
non-MUSCL  upwind  TVD  scheme  of  Yee.^  These  numerical  techniques  were  developed 
for  one-dimensional,  perfect-gas  flows  but  have  been  also  applied  to  a  wide  variety  of 
multidimensional,  real  gas  flows  in  thermochemical  nonequilibrium  with  great  success. 
Only  cursory  descriptions  of  Roe’s  scheme  and  the  Yee  TVD  extension  are  given  here; 
thorough  treatment  of  the  underlying  theory  is  well  documented  in  the  references. 

Most  of  the  modern  upwind  schemes  are  based  on  the  signs  of  the  eigenvalues  of  the 
flux  Jacobians.  The  Jacobian  A  is  decomposed  to 

A  =  (3.4.1) 

where  is  a  diagonal  matrix  with  the  eigenvalues  of  A  on  the  diagonal  and  P^  is  the 
corresponding  matrix  of  right  eigenvectors.  In  practice  it  is  difficult  to  obtain  these  matri¬ 
ces  unless  a  similarity  transform  is  used.  In  that  case, 

A  =  MAM~^  =  MP~A-P~^M~^  (3.4.2) 

AAA 

where 


1.  See  Roe  (1984). 

2.  See  Yee  (1989). 


(3.4.3) 
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^  =  w 


(3.4.4) 


The  eigenvalues  of  A  and  A  are  the  same.  The  vector  V  is  composed  of  the  primitive  vari¬ 


ables: 


(Pl»  •••  >  Pne»  Pne+l’ ’  P/^^’ • 


The  same  decomposition  is  done  for  the  r\  direction, 


B  =  MBM  ^  =  MP-A-P-^M  ^  =  Pj^A^pJ . 

B  B  B  ts  D  IS 


(3.4.5) 


(3.4.6) 


A  detailed  explanation  of  the  diagonalization  process  and  resulting  definition  of  the  above 
matrices  is  given  in  Appendix  D.^ 

Given  the  above  matrices,  many  different  upwind  schemes  can  be  implemented.  Gen¬ 
erally  they  are  differentiated  by  the  choice  of  the  ‘average  state’  at  which  the  Jacobian 
terms  are  evaluated.  The  average  state  for  the  i  +  1/2  cell  face,  denoted  here  by  the  sub¬ 
script  avg,  is  some  weighted  average  of  the  i  and  i  +  1  states.  Roe’s  approximate  Riemann 
solver  represents  the  flux  as 


F' 

^  I  1+1/2 


where 


2  I  ^  I  i+1  ^A,  avgl'^A,  avg  ®A,  avg^ 


=  Plavg(Uu,-U,). 


(3.4.7) 


(3.4.8) 


In  practice  the  components  that  comprise  the  eigenvalue  matrix  are  limited  to  pre¬ 

vent  the  appearance  of  aphysical  expansion  shocks.  The  superscript  (1)  refers  to  the  partic¬ 
ular  components  of  the  vector.  The  are  replaced  by  the  function  \)/(  '^^^avg^ ’ 


where 


3.  Thorough  treatment  of  the  relationships  between  the  conserved,  primitive,  and  characteristic 
variables  for  nonreacting  flows  is  given  in  Hirsch  (1990),  Chapter  16.  Extension  to  reacting 
flows  is  straightforward. 
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¥(z)  = 


(3.4.9) 


kl  ^  8i 
kl  ^  Sj 

and 

6i  =  82 (|a'|  +  |v'|  +  |a  (5./  +  5./)  |)  .  (3.4.10) 

In  the  original  reference,^  suggested  values  for  82  are  between  0.05  and  0.25.  In  this 
work,  82  is  specified  as  a  function  of  the  pressure  gradient: 

82  =  max  (0.15, 0.5  |pg|)  ,  (3.4.11) 

where 
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,  2  ~2 

z  +8j 
.  ^8^ 


P8  = 


Pmj-Pu 


(3.4.12) 


This  specification  of  82  increases  numerical  damping  in  the  vicinity  of  shocks. 

For  Roe’s  scheme,  the  average  state  is  chosen  such  that  the  numerical  flux  exactly  cap¬ 
tures  shock  and  contact  discontinuities.  The  Roe-averaged  variables  are  defined  as  fol¬ 
lows: 


and 


-  ^  Pe,i  +  ^Pe,M 

^e,  avg  1  +  r 
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1  +  r 
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(3.4.13) 

(3.4.14) 

(3.4.15) 


(3.4.16) 


where 


4. 


See  Yee  (1989). 
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The  total  enthalpy  is  defined  by 


jj^2E±2,  (3.4.18) 

P 

All  other  variables  at  the  averaged  state  are  determined  from  the  above  quantities.  In  a 
similar  fashion  Roe’s  scheme  is  applied  to  the  y  +  1/2  cell  face  to  obtain  G' j  y  +  ]/2  • 

Only  slight  modifications  are  required  for  the  TVD  extension  to  second-order  accu¬ 
racy.  The  non-MUSCL,  upwind,  TVD  numerical  flux  of  Yee  is 


^  I  i+l/2  ~  2^^  I  I  '+1  avg^A,  avg^  ’ 


(3.4.19) 


where  the  average  state  can  be  determined  from  any  of  the  first-order  upwind  schemes.  In 
this  work,  the  Roe  average  is  used.  The  components  of  the  vector  O  are  denoted  by  (j) 
and  are  defined  by 


^  i4,  avg  Ta,  avg 

The  function  \\f,  the  eigenvalues  X 
flows,  the  function  a  is  defined  by 


+*/]  ■  (3-4-20) 

,  and  the  are  the  same  as  before.  For  steady 


cy{z)  =  (3.4.21) 

The  is  the  limiter.  The  limiter  controls  how  much  numerical  dissipation  is  added  to 
damp  spurious  oscillations  associated  with  second-order  accuracy.  For  all  fields  (all  (/)),  a 
minmod  limiter  used  in  the  current  work: 


where 


.(  (9  (9 


(3.4.22) 


minmod  {x,  y)  =  sign  (x)  ■  max  (0,  min  [\x\,y  •  sign  (x)])  .  (3.4.23) 

The  subscripts  on  refer  to  the  average  states  for  the  cell  faces  on  the  left  (  / -  1/2) 
and  right  ( i  +  1  /2 )  sides  of  cell  i.  Other  limiters  were  implemented,  but  no  significant  dif¬ 
ferences  were  seen  in  the  solutions.  Finally, 
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(3.4.24) 


This  completes  the  discussion  of  the  explicit  inviscid  terms. 


The  terms  comprising  the  viscous  fluxes  are  centrally  differenced  because  they  are 
parabolic  in  nature.  The  evaluation  of  F'y  1  +  1/2  G'y  j+1/2  requires  property  values 
and  derivatives  at  the  cell  face  of  interest.  Cell  face  property  values  are  found  by  averag¬ 
ing  the  quantities  at  the  adjoining  cells,  i.e.  p;  +  j /2  =  ( P,-  +  P,+i )  .  Two  types  of  deriv¬ 
atives  are  required.  Property  derivatives  are  calculated  in  the  computational  domain.  Mesh 
derivatives  relate  changes  in  the  computational  space  to  changes  in  the  physical  space. 

Expressions  for  and  are  needed,  where  (p  is  any  of  the  properties  ,  u,  v,  or  T. 
Given  these  terms  F'y  1+1/2  ^ V7  +  1/2  calculated  from  Equations  (3.3.7)  and 

(3.3.8).  Applying  the  chain  rule. 


and 
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(3.4.26) 


As  discussed  in  Appendix  C,  the  mesh  derivatives  (or  mesh  metrics)  ||,  and  ^ 

are  more  easily  calculated  in  terms  of  their  inverses.  Then 


and 


where 
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Consider  the  i  +  1  /2  cell  face.  Property  derivatives  in  the  ^  -direction  are  centrally 
differenced  by 


dtp 


lj+l/2,; 


while  those  in  the  r\  -direction  are  written 


(3.4.30) 
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dp 
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For  the  7+1/2  cell  face  the  following  formulae  are  used: 


dtp 


/,7+l/2 
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and 


(3.4.31) 
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(3.4.33) 


With  the  above  expressions,  the  terms  comprising  the  viscous  fluxes  can  be  calculated 
using  the  relations  in  Chapter  II. 


The  relations  introduced  in  section  n.6  are  used  to  calculate  the  net  production  of  the 
retained  species  and,  if  applicable,  the  source  terms  arising  from  the  axisymmetric  coordi¬ 
nate  system.  The  evaluation  of  the  axisymmetric  terms  is  straight-forward.  With  slight 
modifications,  the  chemical  source  terms  are  calculated  using  the  relations  described  in 
Chapter  E. 

For  some  combusting  flowfield  simulations,  it  is  possible  for  a  fast  chemical  reaction 
to  completely  convert  reactants  to  products  within  a  single  mesh  cell;  that  is,  a  reaction 
reaches  equilibrium  in  a  fraction  of  the  time  it  takes  for  the  reacting  molecules  to  cross  the 
cell.  The  energy  released,  and  thus  the  temperature,  is  then  overpredicted.  Since  reaction 
rates  increase  with  temperature,  the  overprediction  is  amplified  on  the  next  iteration.  The 
overall  result  is  that  ‘runaway’  reactions  occur  and  the  reaction  front  moves  at  an  aphysi- 
cal  speed.^  It  is  pointed  out  by  Eberhardt  et.  al.^  that  while  the  grid  cannot  resolve  such 
fast  reactions,  the  associated  difficulties  can  be  circumvented  by  limiting  the  Damkohler 
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number.  This  is  done  by  calculating  a  time  scale  for  each  reaction  and  a  fluid  time  scale, 
comparing  them,  and  then  limiting  the  reaction  rates  of  the  disagreeable  reactions. 

The  Damkohler  number  is  the  ratio  of  the  fluid  time  scale  to  the  time  scale  of  the  reac¬ 
tion,  defined  as 


Da 

X 

rxrii  r 

where  r  is  the  index  of  the  reaction.  The  fluid  time  scale  for  a  cell  is  defined  as 

(\u'\+a)S.+  (|v'hfl)5. 

T  =  _ t _ I 

fluid  Q 


(3.4.34) 


(3.4.35) 


The  reaction  time  scales  have  different  forms  depending  on  the  number  of  reactant  and 
product  species.  Time  scales  for  the  forward  and  backward  components  of  a  reaction  are 
calculated,  and  the  more  limiting  of  the  two  is  used.  For  the  generic  exchange  reaction 

A  +  5  C  +  Z) , 

the  two  time  scales  are 


and 

(3.4.37) 

These  time  scales  are  estimates  of  the  time  required  for  reactant  concentrations  to  decay  to 
1/e  of  the  initial  concentrations.^  For  the  forward  reaction  A  and  B  are  the  reactants,  but 
C  and  D  are  the  reactants  in  the  backward  reaction.  For  other  reaction  forms,  the  time 
scales  are  slightly  modified  but  are  found  in  the  same  manner. 

Once  the  above  time  scales  are  calculated,  the  forward  and  backward  reaction  rate 
constants  are  multiplied  by  the  factor 

fr  =  (3.4.38) 

where 


5.  See  Oran  and  Boris  (1987),  pp.  523-524. 

6.  See  Eberhardt  and  Imlay  (1989). 

7.  See  Kuo  (1986),  p.  154. 
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min(x  f  ,Da  ..) 
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(3.4.39) 
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(3.4.40) 

t  u 
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The  Damkdhler  number  is  limited  to  ,  which  in  this  work  is  10.  For  Da  =  100 , 
no  differences  were  observed  in  the  solutions  but  convergence  was  degraded.  Both  the  for¬ 
ward  and  backward  rate  coefficients  are  multiplied  by  the  same  factor,  so  the  equilibrium 
constant  is  not  affected  by  the  limiting. 

1II.5  Evaluation  of  the  Implicit  Terms 

The  methods  for  calculating  the  explicit  fluxes  F'"  and  G'”  have  been  described. 
Solution  advancement  schemes  have  been  developed  that  require  only  these  explicit  fluxes 
and  are  aptly  called  explicit  schemes.  The  time  step  for  explicit  schemes  is  limited  by  sta¬ 
bility  considerations,  and  for  flows  with  fast  chemical  reactions  the  limitation  can  be 
severe.  Implicit  schemes  alleviate  the  strictures  on  the  time  step,  but  the  cost  per  iteration 
increases.  The  additional  expense  is  for  the  calculation,  and  in  particular  the  inversion,  of 
the  Jacobians. 

The  terms  on  the  explicit  side  are  modelled  as  accurately  as  is  reasonable,  because  the 
final  numerical  solution  satisfies  the  physics  described  by  the  right-hand  side.  In  contrast, 

H  ft  M 

the  implicit  terms,  namely  A  ,B  ,  and  Z  ,  are  usually  approximated.  Previously  the  flux 
Jacobians  were  decomposed  for  the  implementation  of  Roe’s  scheme,  but  they  are  evalu¬ 
ated  differently  for  the  implicit  side.  If  properly  designed  and  implemented,  an  implicit 
scheme  will  do  two  things.  First,  it  will  converge  to  the  numerical  solution  dictated  by  the 
right-hand  side.  Second,  it  will  reduce  the  number  of  iterations  required,  with  the  ultimate 
goal  being  the  reduction  of  overall  computer  time.  The  approximations  made  on  the  left- 
hand  side  distinguish  the  different  implicit  solution  advancement  methods.  More  accurate 
approximations  generally  reduce  the  number  of  necessary  iterations,  but  the  cost  per  time 
step  increases.  The  best  solution  advancement  scheme,  in  terms  of  overall  computer  time. 
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depends  on  the  particular  flow  problem  and  machine  architecture. 

A  modifled  form  of  the  diagonal  implicit  scheme  of  Eberhardt  and  Imlay^  is  used  in 
this  thesis.  The  diagonal  implicit  scheme  is  an  extension  of  the  Lower-Upper-Symmetric- 
Gauss-Seidel  (abbreviated  LU-SGS)  method  of  Yoon  and  Jameson^  to  chemically  reacting 
flows.  To  apply  the  LU-SGS  scheme  the  fluxes  on  the  implicit  side  are  upwinded: 


/  is,  ,A" 


(3.5.1) 


111  1 


The  variables  y.  and  are  averages  of  the  values  at  the  « +  1/2  and  /  -  1/2  cell  faces. 
Consequently  y .  =  y^  .  . .  The  flux  Jacobians  are  discussed  in  Appendix  D.  The  split  flux 
Jacobians  are  constructed  such  that  A”  contains  strictly  nonnegative  eigenvalues  and  a" 
contains  strictly  nonpositive  eigenvalues.  Only  the  inviscid  terms  are  included  in  the  Jaco¬ 
bians.  In  a  similar  fashion  is  split.  As  a  result  Equation  (3.3.15)  is  rewritten 

^~i 
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(3.5.2) 


.  ,  =  -AC/  . 

;-i  hj 


The  LU-SGS  scheme  represents  the  split  flux  Jacobians  as  follows: 

Aj=i(A"±v]  B"  =  i(B“±r5/].  (3.5.3) 

The  spectral  radii  of  A”  and  b”  are  and  ,  respectively.  With  this  representation  of 
the  split  flux  Jacobians,  note  that 


and 


(3.5.4) 
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Consequently  Equation  (3.5.2)  can  be  simplified  to 


(3.5.5) 


8.  See  Eberhardt  and  Imlay  (1989). 

9.  See  Yoon  and  Jameson  (1987)  and  (1988). 


40 


/  i  i^i  +  i  i 
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(3.5.6) 
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Observe  that  for  the  nonreacting  case  the  chemical  Jacobian  is  zero,  and  the  terms  in 
square  brackets  form  a  diagonal  matrix.  The  inversion  of  a  diagonal  matrix  is  much  less 
expensive,  computationally,  than  a  full  matrix,  and  for  this  reason  the  LU-SGS  scheme  has 
a  very  low  cost-per-time-step  compared  to  other  implicit  methods. 

The  spectral  radii  are  approximated  by  the  maximum  eigenvalues: 


=  |m  +fl| 


r^  =  \v  +a\. 


(3.5.7) 


Although  the  viscous  terms  are  not  included  in  the  split  flux  Jacobians,  an  adjustment  is 
made  to  the  spectral  radii  for  viscous  flows.  The  term 

^ - 1  (3.5.8) 


is  added  to  ,  where  C^is  defined  in  Appendix  D.  The  same  term  is  added  to  except 
S-  replaces  5  • . 

With  the  inclusion  of  the  full  chemical  Jacobian,  the  off-diagonal  terms  are  no  longer 
zero  and  the  full  matrix  must  be  inverted.  Eberhardt  and  Imlay  approximated  z”^.  in  such 
a  way  that  the  diagonal  form  is  retained.  In  the  diagonal  implicit  scheme,  only  the  first  ns 
rows  and  ns  columns  of  the  chemical  Jacobian  are  considered.  The  matrix  z”  .  is  repre¬ 


sented  by 


z"^.  =  diagri  0,0, o'). 


(3.5.9) 


Each  chemical  species  has  a  chemical  time  scale,  ,  that  is  based  on  the  reactions  in 
which  the  species  participates.  The  time  scale  is  defined  by 


1  JvtelT 

T  “  Pc  2w 


(3.5.10) 
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where  ^  are  the  components  of  the  first  ns  rows  and  columns  of  the  full  chemical  Jaco¬ 
bian. 

Based  on  its  corresponding  time  scale,  each  species  equation  evolves  at  a  different 
rate.  The  species  equations  are  underrelaxed  relative  to  the  momentum  and  energy  equa¬ 
tions.  The  parameter  [3^  is  an  additional  underrelaxation  factor.  For  a  species  that  partici¬ 
pates  in  fast  chemical  reactions,  1/t^  is  large;  the  corresponding  mass  conservation 
equation  is  underrelaxed  to  a  greater  degree  than  the  equations  associated  with  smaller 
values  of  l/x^ .  For  most  other  implicit  methods,  all  equations  evolve  at  the  same  rate.  An 
advantage  of  the  diagonal  implicit  scheme  is  that  all  equations  are  not  slowed  to  a  time 
step  limited  by  the  most  restrictive  equation;  each  equation  evolves  at  a  rate  determined 
independently  of  the  other  equations. 

The  disadvantages  of  the  scheme  are  also  a  result  of  the  differing  relaxation  rates.  As 
detailed  by  Candler  and  Olynick^®  and  Hassan  et.  al.,^^  the  convergence  properties  of  the 
diagonal  implicit  scheme  suffer  from  several  inconsistencies  caused  by  the  treatment  of 
the  chemical  Jacobian.  There  are  three  disadvantages.  First,  elemental  mass  is  not  con¬ 
served  when  species  evolve  at  different  rates.  Second,  since  the  total  density  is  found  by 
summing  the  species  densities,  it  evolves  at  a  rate  that  is  a  weighted  sum  of  the  species 
evolution  rates.  Compared  to  the  momentum  and  energy  equations,  the  species  equations 
are  underrelaxed.  The  total  density  therefore  lags  the  momentum  and  energy.  The  above 
authors  have  mitigated  these  difficulties  by  replacing  some  of  the  species  mass  conserva¬ 
tion  equations  by  elemental  mass  conservation  equations.  The  chemical  Jacobian  of  this 
modified  diagonal  implicit  scheme  is  represented  by 


Z”.  =  diagfo,...,0,-^,...,;;^-,0,0,ol 

V  ^ne+l  ^ns  ^ 
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10.  See  Candler  and  Olynick  (1991). 

11.  See  Hassan,  Candler,  and  Olynick  (1992). 
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The  underrelaxation  parameter  was  found  to  be  unnecessary  using  the  elemental  for¬ 
mulation.  The  elemental  equations  have  l/x^  =  0  because  there  are  no  source  terms  for 
the  elemental  densities.  When  the  total  density  is  determined  by  summing  the  elemental 
densities,  it  consistent  with  the  momentum  and  energy  because  the  elemental  mass, 
momentum,  and  energy  equations  are  all  releixed  at  the  same  rate.  Because  of  these 
improvements,  the  elemental  conservation  equations  are  employed  in  this  work. 

The  third  disadvantage  is  that  the  diagonal  implicit  scheme,  in  both  the  original  and 
modified  forms,  is  suitable  for  finding  steady-state  solutions  only.  Temporal  accuracy  is 
forfeited  because  each  species  equation  effectively  has  a  different  time  step.  For  each  spe¬ 
cies  equation  there  is  a  different  transformation  between  physical  time  and  computational 
time  (iteration  number).  Thus,  there  is  no  way  to  relate  a  single  physical  time  to  a  given 
iteration  number. 

While  the  elemental  formulation  relieves  some  of  the  difficulties  of  the  diagonal 
implicit  scheme,  another  is  created.  Specifically,  it  is  possible  to  calculate  negative  densi¬ 
ties  for  the  replaced  species.  The  replaced  species  are  those  omitted  from  U,  the  vector  of 
conserved  variables;  they  have  been  replaced  by  elemental  densities.  The  species  included 
in  the  conserved  variable  vector  are  referred  to  as  the  retained  species.  The  replaced  spe¬ 
cies  densities  are  computed  from  the  elemental  and  retained  species  densities  using  Equa¬ 
tion  (2.2.21).  That  equation  is  valid  for  a  specific  location  in  physical  space  (jc,  y)  and 
time  (0  .  However,  it  is  implemented  for  each  point  in  computational  space  (^,1))  and 
time  (iteration  number  n)  .  As  mentioned  previously,  at  a  given  iteration  the  retained 
species  are  at  different  points  in  physical  time  because  each  retained  species  equation  is 
relaxed  at  a  different  rate.  Strictly  speaking,  the  linear  combination  for  computing  the 
replaced  species  does  not  apply.  As  an  simplistic  example  of  how  negative  densities  arise, 
suppose  the  species  CO  and  CO2  are  retained  while  CH^  is  replaced.  Further  assume 
that  CO  has  a  small  chemical  time  scale  compared  to  that  of  CO^  •  If  CO  is  being 
oxidized  to  CO^ ,  the  production  of  CO^  can  be  far  greater  than  the  depletion  of  CO 
because  the  latter  has  a  smaller  time  scale.  When  elemental  conservation  is  enforced  (for 
carbon,  in  this  example),  the  C  in  CO^  must  come  from  the  replaced  species  CH^ .  A 
negative  density  is  computed  when  more  C  is  taken  from  CH^  than  actually  exists  even  if 
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that  is  not  physically  possible. 

To  this  point  the  equations  have  been  developed  without  specifying  which  of  the  spe¬ 
cies  are  retained  and  which  are  replaced.  The  development  of  the  flowfield  from  the  initial 
conditions  to  the  final  numerical  solution  depends  on  the  choice  of  retained  species, 
because  only  the  ne  +\  to  ns  time  scales  are  used.  For  some  selections  negative  densities 
may  arise,  while  for  others  all  replaced  densities  will  be  positive.  Based  on  careful  numer¬ 
ical  experiments,  the  replaced  species  are  CH^ ,  O^,  H2,  and  N2  in  this  thesis.  No  signif¬ 
icant  negative  densities  were  observed  in  the  flowfields  presented. 

The  occurrence  of  negative  densities  was  not  reported  by  Candler  and  Olynick  or  by 
Hassan  et.  al.,  who  used  the  elemental  formulation  for  reentry  flowfields.  The  gas  model  in 
those  studies  was  8-  or  11 -species  air,  both  of  which  consider  only  the  elements  nitrogen 
and  oxygen.  It  is  possible  that  negative  densities  are  more  likely  to  occur  for  gas  models 
with  more  elements.  It  is  also  possible  that  they  arise  because  of  the  global  reaction,  which 
has  a  nonstandard  rate  form.  However,  the  example  above  suggests  chemical  models  with 
only  Arrhenius  rates  are  no  less  susceptible. 

The  diagonal  implicit  scheme  is  not  time-accurate,  and  it  could  be  argued  that  there¬ 
fore  only  the  steady  state  solution  has  any  physical  meaning.  The  occurrence  of  negative 
densities  during  the  transient  should  not  cause  alarm  because  the  intermediate  solutions 
have  no  significance.  However,  the  use  of  these  negative  densities  in  many  of  the  relations 
described  in  Chapter  II  can  lead  to  inconsistencies  that  ultimately  prevent  the  development 
of  the  correct  numerical  solution.  During  testing  of  the  algorithm,  some  choices  for 
replaced  species  resulted  in  negative,  monotonically  decreasing  CO 2  mass  fractions. 
Obviously  the  final  solution  was  aphysical.  For  this  reason,  thorough  testing  and  valida¬ 
tion  of  computer  programs  employing  the  diagonal  implicit  scheme  is  highly  recom¬ 
mended.  Clearly  further  investigation  into  the  conditions  which  lead  to  negative  replaced 
species  densities  is  warranted. 

1II.6  Solution  Advancement 

When  the  entire  computational  domain  is  considered.  Equation  (3.5.6)  is  a  large, 
block-pentadiagonal  system.  At  the  beginning  of  an  iteration,  an  approximate  solution  for 


the  entire  flowfield  is  known.  The  solution  is  improved  during  the  iteration,  and  after  a 
number  of  iterations  it  is  a  close  approximation  to  the  physical  solution  described  by  the 
governing  equations. 

An  iteration  consists  of  several  steps.  First  the  explicit  terms  are  evaluated,  as  dis¬ 
cussed  in  Section  III.4,  and  Af/”  is  calculated  for  the  interior  mesh  cells.  The  interior  cells 
are  all  those  within  the  boundary  of  the  physical  domain.  For  a  grid  of  size  il  by  jl,  where 
il  is  the  number  of  points  in  the  ^  -direction  and  jl  is  the  number  in  the  r\  -direction,  the 
interior  cells  have  i  =  2, ...,  il-  1  and  j  =  2,  ...Jl-  1 . 

Second,  Equation  (3.5.6)  is  approximately  solved  by  two  Gauss-Seidel  sweeps.  On  the 
forward  sweep,  from  the  lower  left  corner  to  the  upper  right  corner  of  the  mesh, 

i  1  i^i  iK.  i  1  /*  •  (3.6.1) 

lyj  «-l  +,  l-l,J  ^  ' 
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On  the  backward  sweep, 

+  (3.6.2) 

For  both  sweeps,  the  diagonal  matrix  is 


D.  ■  =  y  .  . 
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it  is  this  diagonal  matrix  that  is  easily  inverted,  resulting  in  a  low  cost  per  iteration. 

At  the  end  of  the  backward  sweep  6c/*  is  known  for  the  interior  cells.  The  third  step  is 
the  update  step: 


u"*'  =  v"  +  iu".  (3.6.4) 

Any  other  necessary  flowfield  properties  are  determined  from  the  updated  conserved  vari¬ 
ables.  The  boundary  cells  are  updated  as  the  fourth  step,  which  is  discussed  in  the  next 
section. 


111.7  Boundary  Conditions 

The  cell  face  shared  by  a  boundary  cell  and  the  first  interior  cell  is  the  physical  bound¬ 
ary.  It  may  be  a  solid  surface,  or  it  may  be  a  line  in  space  defining  the  edge  of  the  flowfield 


to  be  simulated.  The  boundary  cells  provide  a  means  of  imposing  boundary  conditions  at 
the  physical  boundaries  and  are  defined  such  that  the  physical  boundary  is  halfway 
between  the  cell  centroids  of  the  boundary  cell  and  first  interior  cell.  Additionally  the 
boundary  and  first  interior  cells  have  the  same  volumes. 

At  an  inflow  boundary,  the  conditions  are  specified  by  the  freestream  conditions.  At  an 
outflow  boundary  zeroth-order  extrapolation  is  applied: 


%,.J  =  (3-7-1) 

where  the  boundary  cell  has  indices  il,  j,  and  (p  is  any  property.  Although  this  extrapola¬ 
tion  is  strictly  valid  for  supersonic  outflow  only,  it  is  also  applied  where  boundary  layers 
cross  the  outflow  boundary. 

For  inviscid  simulations,  the  flow  tangency  condition  specifies  the  boundary  cell  prop¬ 
erties  at  solid  surfaces.  The  same  conditions  are  applied  at  centerline  boundaries.  Consider 
a  boundary  cell  with  indices  l,j.  The  u'  velocity  is  parallel  to  the  boundary,  while  v'  is 
perpendicular  to  it.  The  velocities  at  the  interface  should  be 


u 


2,j 


(3.7.2) 


Solving  for  the  boundary  cell  velocities. 


and 


(3.7.3) 


From  the  definition  of  contravariant  velocities,  «,  ■  and  v,  can  be  determined.  It  is 

^9  J  ^9  J 

assumed  that  the  grid  is  orthogonal  near  the  boundary,  and  the  remaining  properties  are 
reflected: 


<Pl  .  =  (pj  ..  (3.7.5) 

where  cp  is  any  of  the  elemental  or  species  densities,  the  temperature,  or  the  pressure.  The 
remaining  variables  are  calculated  from  these. 

For  viscous  simulations,  the  no-slip  condition  is  appropriate  at  solid  boundaries.  In 
that  case 
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SO 


u\  ■  =  1u'  ,,—u' 

l,j  wall  “  2,j 


and 


wall  ^  2,j- 


(3.7.7) 

(3.7.8) 


For  most  cases  the  walls  do  not  move  in  relation  to  the  frame  of  reference  of  the  simula¬ 
tion,  and  and  are  both  zero.  However,  the  accelerator  tube  wall  moves  with 
respect  to  the  projectile  surface,  so  for  viscous  ram-accelerator  simulations 
and  =  0 .  The  remaining  variables  are  reflected,  as  in  the  inviscid  surface  case.  For 
such  a  specification  an  adiabatic  boundary  is  assumed  since  T.  .  =  .. 


47 


CHAPTER  IV 


BLUNT  BODY  SIMULATIONS 


rv.l  Introduction 

The  purpose  of  the  blunt  body  simulations  is  to  validate  the  numerical  method.  The 
simulations  concern  the  hypersonic,  inviscid,  reacting  flow  over  cylinders  of  various 
diameters.  Blunt  body  flowfields  contain  a  bow  shock  that  is  normal  near  the  stagnation 
line  but  weakens  away  from  it.  A  moderate  expansion  is  also  present  as  flow  accelerates 
around  the  body.  Accurately  capturing  these  features  requires  a  robust  flow  solver. 

For  combusting  flows,  the  bow  shock  prompts  ignition  because  of  the  high  tempera¬ 
tures  in  the  shock  layer.  Ignition  generally  refers  to  the  beginning  of  rapid  energy  release. 
Depending  on  the  conditions,  a  steady  energy-release  front  can  form  between  the  bow 
shock  and  the  body.  The  energy-release  front  is  also  called  the  reaction  front  because  it  is 
where  the  fuel  is  consumed  and  converted  to  combustion  products.  The  region  between 
this  front  and  the  shockwave  is  called  the  induction  zone;  after  reaching  a  temperature 
high  enough  to  initiate  ignition,  there  is  a  time  delay  before  the  energy  is  released,  called 
the  induction  time.  The  length  of  the  induction  zone  is  essentially  the  distance  covered  by 
the  reacting  particles  during  the  induction  time. 

Accurate  modelling  of  the  induction  zone  requires  a  sound  chemical  kinetics  mecha¬ 
nism.  The  induction  time  is  highly  dependent  on  the  temperature,  pressure,  and  composi¬ 
tion  of  the  gas.  Thus,  the  shock  and  reaction  front  standoff  distances  are  sensitive  to  the 
chemistry  model  implemented,  and  for  this  reason  are  used  as  the  basis  of  comparison  to 
established  solutions.  The  results  exhibit  the  advantages  and  drawbacks  of  the  quasi-glo- 
bal  mechanism  and  demonstrate  the  ability  of  the  numerical  scheme  to  simulate  hyper¬ 
sonic,  combusting  flowfields. 

rv.2  Flow  Conditions 

The  blunt  body  flowfields  presented  model  the  flow  over  rods  of  1,  3,  and  7  mm  diam¬ 
eter.  The  gas  mixture  is  stochiometric  methane-air  at  freestream  temperature 
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=  295  K,  pressure  =  51000  Pa,  and  velocity  =  2330  m/s.  The  resulting 
Mach  number  is  6.61.  The  induction  time  depends  only  on  the  conditions  of  the  gas 
behind  the  shockwave.  However,  since  the  bow  shock  and  energy-release  front  standoff 
distances  are  coupled,  the  flowfields  induced  by  each  rod  size  are  not  the  same  when  non- 
dimensionalized  by  diameter. 

The  computational  grid  has  the  same  proportions  for  the  different  rod  sizes;  it  is  lin¬ 
early  scaled  to  the  appropriate  diameter.  To  reduce  the  expense  of  the  simulations,  only 
one  half  of  the  rod  is  considered.  The  full  flowfield  can  be  generated  by  reflection  across 
the  centerline.  Since  the  data  of  interest  is  contained  in  the  flow  ahead  of  the  body,  the 
wake  region  is  not  included  in  the  simulation.  The  mesh  is  shown  in  Figure  4.1.  The  grid 
consists  of  95  points  radially  and  91  points  around  the  body.  The  boundary  conditions  are 
freestream  conditions  along  the  left  (curved)  boundary,  centerline  conditions  along  the 
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Figure  4.1  Blunt  body  mesh  and  representative  temperature  contours. 


bottom  boundary,  inviscid  wall  along  the  body  surface,  and  outflow  conditions  at  the  right 
edge  of  the  domain.  Additionally  a  temperature  contour  plot  is  shown  in  Figure  4.1.  The 
temperature  contours  illustrate  the  relevant  flowfield  features.  The  shockwave  and  energy- 
release  front  are  unambiguous. 

The  cases  were  originally  investigated  by  Yungster  and  Rabinowitz.^  Those  authors 
used  a  20  species,  52  reaction  chemical  kinetics  mechanism,  which  was  carefully  derived 
from  a  33  species,  149  reaction  model^  based  on  a  comparisons  with  experimentally  mea¬ 
sured  ignition  delay  times  in  shock  tubes.  Because  of  its  systematic  development,  the 
Yungster  and  Rabinowitz  solutions  are  the  benchmark  for  the  validations  in  the  present 
research.  The  data  for  comparison  are  pressure  and  temperature  profiles  on  the  centerline. 
Pressure  profiles  clearly  illustrate  the  shock  standoff  distance,  while  the  temperature  pro¬ 
files  additionally  show  the  energy-release  front. 

IV.3  Discussion  of  the  Simulations 

Flowfield  simulations  were  obtained  for  the  conditions  described  above.  Two  sets  of 
runs  were  made.  In  the  first  set,  the  value  of  the  pre-exponential  coefficient  Ag  is  that  sug¬ 
gested  in  the  original  reference.^  In  the  second  set,  double  the  original  value  is  used. 
Recall  that  the  subscript  6  refers  to  the  sixth  reaction,  which  is  the  global  reaction.  As  dis¬ 
cussed  in  section  II.6,  the  pre-exponential  coefficient  should  be  calibrated  to  match  estab¬ 
lished  data. 

Each  reacting  simulation  ran  between  4000  and  6000  iterations.  The  nonreacting  run 
took  2000  iterations.  For  all  cases  the  residual  dropped  at  least  4  orders  of  magnitude.  In 
some  instances  an  interaction  between  the  Damkohler  limiting  and  the  global  reaction  rate 
inhibited  further  reduction  of  the  residual.  However,  flowfield  properties  remained 
unchanged  after  ~  2500  iterations. 

For  the  results  obtained  using  the  original  rate,  the  nondimensionalized  pressure  pro¬ 
files  along  the  centerline  are  shown  in  Figure  4.2.  The  pressure  is  nondimensionalized  by 


1 .  See  Yungster  and  Rabinowitz  (1993). 

2.  See  Frenklach,  Wang,  and  Rabinowitz  (1992) 

3.  See  Westbrook  and  Dryer  (1981). 


and  the  x-axis  is  nondimensionalized  by  the  radius  of  the  rod.  The  established  data  of 
Yungster  is  represented  by  dashed  lines,  while  for  the  present  results  they  are  solid.  The 
pressure  profile  for  the  nonreacting  case  is  also  shown  for  comparison  and  is  denoted  by 
the  solid  line  with  square  symbols.  Flow  is  from  left  to  right.  For  all  cases  the  bow  shock 
is  indicated  by  the  steep  pressure  rise.  It  is  clear  from  the  plot  that  the  shock  standoff  dis¬ 
tance  is  underpredicted  using  the  original  rate  coefficient  . 

The  nondimensional  temperature  profiles  are  shown  in  Figure  4.3  with  the  same  line 
and  symbol  conventions.  The  first  increase  in  temperature  denotes  the  shockwave.  The 
second  increase  signifies  the  energy-release  front.  The  benchmark  profiles  show  the  reac¬ 
tion  front  further  from  the  body  than  the  current  results.  On  the  basis  of  these  compari¬ 
sons,  the  conclusion  is  that  the  original  value  of  the  pre-exponential  coefficient  is  too  low. 
If  the  induction  time  was  smaller,  the  standoff  distance  of  the  energy-release  front  would 
be  further  from  the  body.  As  a  second  consequence,  the  shock  standoff  distance  would 
also  be  greater. 

A  second  set  of  runs  was  generated  with  double  the  original  value  for  .  The  corre¬ 
sponding  pressure  and  temperature  profiles  are  shown  in  Figures  4.4  and  4.5,  respectively. 
Clearly  the  pressure  and  temperature  profiles  are  in  better  agreement  with  the  doubled 
rate.  The  induction  zone  is  roughly  the  same  length  for  the  present  results  and  the  simula¬ 
tions  of  Yungster  and  Rabinowitz,  indicating  the  different  kinetics  models  are  in  agree¬ 
ment. 

Although  the  shock  and  energy-release  front  standoff  distances  and  the  induction  zone 
length  are  improved  using  the  doubled  rate,  there  are  still  slight  incongruities  in  the 
results.  While  the  profiles  are  similar,  their  shapes  are  not  the  same.  Compared  to  the  base¬ 
line  data,  the  energy-release  front  is  smeared.  The  temperature  behind  the  energy-release 
front  is  higher.  The  shock  standoff  distance  matches  very  well  for  the  1  mm  rod,  but  dis¬ 
crepancies  can  be  seen  for  the  larger  diameters.  The  results  underscore  the  shortcomings 
of  the  quasi-global  mechanism;  by  modelling  many  detailed  reactions  with  a  single  phe¬ 
nomenological  reaction,  the  details  of  the  reaction  front  are  sacrificed. 

On  the  other  hand,  the  advantage  of  the  quasi-global  model  is  that  reasonable  esti¬ 
mates  of  the  induction  zone,  standoff  distances,  and  gas  composition  and  temperature  are 
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Figure  4.3  Temperature  profile  comparison  using  the  original  rate. 


Figure  4.5  Temperature  profile  comparison  using  the  doubled  rate. 


obtained  at  a  lower  computational  cost.  For  the  quasi-global  model,  13  partial  differential 
equations  are  solved  at  each  grid  point:  10  elemental  or  species  mass  equations,  2  momen¬ 
tum  equations,  and  1  energy  equation.  If  the  mechanism  of  Yungster  and  Rabinowitz  were 
implemented,  there  would  be  10  more  mass  conservation  equations,  bringing  the  total  to 
23.  The  computational  cost  per  iteration  scales  roughly  as  the  square  of  the  number  of 
equations.  Thus,  doubling  the  number  of  species  almost  quadruples  the  expense  of  the 
solution. 

The  outcome  of  the  blunt  body  simulations  is  that  the  numerical  method  is  validated. 
Conditions  are  representative  of  those  in  ram-accelerator  flows  in  that  strong  shock  waves 
and  expansions  are  present,  and  the  freestream  Mach  number,  temperature,  pressure,  and 
mixture  composition  are  similar  to  those  of  the  simulations  presented  in  the  next  chapter. 
With  the  doubling  of  the  pre-exponential  coefficient  of  the  global  reaction,  the  induction 
time  and  associated  flowfield  effects  are  accurately  approximated.  While  more  detailed 
reaction  mechanisms  exist,  the  quasi-global  mechanism  is  suitable  at  a  considerably  lower 
computational  cost. 


CHAPTER  V 


RAM- ACCELERATOR  SIMULATIONS 


V.l  Introduction 

In  this  chapter  the  ram-accelerator  flowfields  are  presented.  The  purpose  of  these  sim¬ 
ulations  is  to  investigate  the  performance  of  the  double-cone  projectile  geometry.  In  par¬ 
ticular,  the  unstart  characteristics  of  such  a  projectile  are  sought.  Stable  ram-accelerator 
operation  may  be  inhibited  for  a  variety  of  reasons,  some  of  which  were  mentioned  in 
Chapter  I.  In  this  thesis  only  gasdynamic  circumstances  are  considered,  and  in  particular 
the  focus  is  on  unstarts  due  to  combustion  on  the  forebody. 

V.2  Ram-Accelerator  Geometry  and  Simulation  Conditions 

For  all  the  simulations  the  solution  technique  described  in  Chapter  III  is  utilized.  The 
flowfields  are  assumed  to  be  axisymmetric  and  laminar  and  the  gas  thermally  perfect.  As  a 
result  of  the  blunt  body  solutions  discussed  in  Chapter  IV,  the  pre-exponential  coefficient 
for  the  global  reaction  is  doubled. 

The  projectile  geometry  of  interest  was  first  reported  by  Yungster  and  Rabinowitz.^ 
The  dimensions  (in  mm)  and  computational  grid  are  shown  in  Figure  5.1.  The  grid  has 
205  X  72  points.  They  are  clustered  in  the  region  near  the  projectile  shoulder  to  resolve  the 
complex  shockwave  interactions  and  chemistry  and  near  the  projectile  and  accelerator 
tube  surfaces  to  resolve  boundary  layers.  Without  adequate  grid  resolution,  the  flowfield 
can  be  substantially  mispredicted.  In  inviscid  simulations  at  Mach  1 1  (not  included  here,) 
a  coarse  mesh  in  the  shoulder  region  resulted  in  a  slightly  weaker  transmitted  shock;  con¬ 
sequently  ignition  did  not  occur  until  behind  the  reflected  shockwave.  On  an  improved 
grid  the  ignition  moved  upstream  to  the  transmitted  shock,  in  agreement  with  the  viscous 
results  presented. 

A  no-slip  condition  is  applied  at  the  projectile  surface  and  the  accelerator  tube  wall. 
Note  that  the  tube  wall  moves  at  a  velocity  with  respect  to  the  projectile.  There  are 


1 .  See  Yungster  and  Rabinowitz  (1993). 
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Accelerator  Tube  Wall 

Figure  5.1  Ram-accelerator  computational  mesh  and  projectile  geometry. 

several  cells  ahead  of  the  projectile  where  a  centerline  boundary  condition  is  imple¬ 
mented.  The  inflow  is  supersonic  and  properties  along  the  inflow  boundary  are  set  to 
freestream  conditions.  Zeroth-order  extrapolation  is  applied  at  the  outflow  boundary. 

The  simulations  do  not  consider  the  region  behind  the  projectile,  referred  to  as  the 
wake,  for  several  reasons.  First  of  all,  extending  the  computational  domain  necessitates 
many  more  grid  points.  In  some  region  behind  the  projectile  base,  flow  is  subsonic.  For 
the  lower  projectile  speeds,  the  subsonic  region  can  extend  many  projectile  lengths  down¬ 
stream.  The  entire  subsonic  region  must  be  included  within  the  domain  for  the  combined 
system  of  partial  differential  equations  and  boundary  conditions  to  be  well-posed.  Such  a 
large  number  of  additional  grid  points  markedly  increases  the  cost  of  the  simulation.  Addi¬ 
tionally  this  subsonic  region  converges  much  more  slowly  than  the  rest  of  the  flowfield 
because  of  its  elliptic  character.  As  a  result,  many  more  iterations  are  required,  also 
increasing  the  cost  of  the  solutions.  Finally,  while  some  mechanisms  of  unstart  rely  on 
thermal  choking  of  the  subsonic  flow  in  the  projectile  wake,  the  current  study  focuses  on 
forebody  preignition  of  the  fuel.  While  the  inclusion  of  the  wake  in  the  simulation  can  be 
valuable,  it  is  not  necessary  for  the  present  investigation. 

The  freestream  gas  is  a  methane-air  mixture  with  equivalence  ratio  of  0.5;  that  is,  there 
is  one-half  of  the  fuel  needed  for  a  stochiometric  mixture.  The  composition  is 
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CH^  +  4  O2  +  15.04  A^2  »  where  nitrogen  is  a  diluent.  The  freestream  conditions  specified 
are  =  1  atm  and  =  300  K.  Five  different  projectile  speeds  are  considered: 
2459.2,  3161.8,  3864.4,  4567.0,  and  5269.6  m/s,  corresponding  to  Mach  numbers  of  7,  9, 
11, 13,  and  15. 

For  the  solutions  presented,  the  residual  is  reduced  at  least  five  orders  of  magnitude. 
Each  simulation  required  between  12,500  and  15,000  iterations,  except  for  the  Mach  9 
case  which  required  22,500  iterations.  Runs  were  carried  out  on  two  machines.  On  a  sin¬ 
gle-processor  Sun  SparcStation  10,  each  1000  iterations  took  11  hours.  The  same  calcula¬ 
tions  require  about  6  hours  on  a  single-processor  SGI  R4400  Indigo. 

V.3  Discussion  of  the  Flowfields 

To  clarify  the  discussion,  the  features  of  the  ram-accelerator  flowfield  are  introduced. 
A  schematic  of  the  flowfield  is  shown  in  Figure  5.2.  The  flow  is  from  left  to  right.  The 
shock  S 1  is  formed  by  the  projectile  tip  and  is  called  the  initial  shock.  The  ramp  shock  is 
labelled  S2  and  results  from  the  increase  in  cone  half-angle.  These  two  shocks  coalesce 
into  the  transmitted  shock,  S3;  after  intersecting  with  the  accelerator  tube  wall,  it  is  called 
the  reflected  shock,  S4.  Behind  the  projectile  shoulder,  flow  is  separated  due  to  the  interac¬ 
tion  between  the  reflected  shock  and  the  boundary  layer  along  the  projectile  surface.  Thor¬ 
ough  discussion  of  this  separation  bubble  is  postponed  to  Section  V.4.  Behind  this 
separation,  flow  reattaches;  as  a  result,  a  reattachment  shock,  S5,  is  formed. 


Figure  5.2  Features  of  the  ram-accelerator  flowfield. 
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The  schematic  illuminates  only  the  shock  system.  Not  shown  are  a  variety  of  contact 
surfaces,  expansions,  energy-release  fronts,  and  their  interactions  with  the  shockwaves 
and  each  other.  For  the  runs  at  lower  projectile  velocities,  there  is  also  a  small  separation 
at  the  beginning  of  the  second  ramp.  Because  of  these  interacting  waves,  the  chemical 
reactions  occurring  in  the  mixture,  and  the  viscous  effects  of  the  boundary  layers  and  sep¬ 
arated  regions,  the  flowfields  presented  are  highly  complex. 

Some  of  the  features  described  are  seen  in  Figure  5.3,  showing  the  Mach  number  con¬ 
tours  for  each  projectile  velocity,  and  Figure  5.4,  showing  the  density  contours.  The  shock 
system  is  clearly  delineated  in  both  figures,  but  the  separation  behind  the  projectile  shoul¬ 
der  is  more  evident  in  the  Mach  number  contours.  The  energy-release  fronts  are  not  appar¬ 
ent;  they  are  revealed  in  mass  fraction  contour  plots.  The  boundary  layer  along  the 
projectile  surface  is  recognizable;  however,  the  accelerator  tube  boundary  layer  is 
obscured  because  the  velocity  gradient  at  the  accelerator  tube  wall  is  much  lower. 


Figure  5.3  Mach  number  contours. 
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For  the  =  7  case,  a  smaller  separation  at  the  beginning  of  the  second  ramp  effec¬ 
tively  rounds  the  compression  corner;  in  place  of  the  S2  shock,  weaker  compression 
waves  form  off  the  separation.  As  the  speed  increases,  the  comer  separation  diminishes  in 
size  until  flow  is  completely  attached  at  =  15 . 

At  the  point  where  the  initial  and  ramp  shocks  combine  to  form  the  transmitted  shock, 
a  contact  surface  divides  the  downstream  flow  into  two  parts.  The  gas  that  has  passed 
through  the  S 1  and  S2  shockwaves  is  not  compressed  as  much  as  that  which  has  passed 
through  the  transmitted  shock.  Thus,  the  contact  surface  divides  gases  at  two  different 
states.  This  state  difference  is  propagated  downstream  and  can  affect  the  ignition  location. 

In  the  projectile  boundary  layer  the  gas  is  heated  significantly.  The  methane  mass  frac¬ 
tion  is  plotted  in  Figure  5.5  and  shows  that  the  boundary  layer  is  hot  enough  to  ignite  the 
fuel  at  all  the  velocities  considered.  For  the  =  7  case,  burning  is  limited  to  the 
boundary  layer.  The  burned  gas  flows  into  the  separation  bubble  and  out  of  the  domain 


M„  =  7 


M„  =  9 

M,  =  11 


M_  =  13 


M,  =  15 


Figure  5.4  Density  contours. 
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without  igniting  the  rest  of  the  flowfield.  Although  gas  is  burning  ahead  of  the  projectile, 
not  enough  energy  is  released  to  cause  unstart. 

For  the  Mach  9  case,  the  reflected  shock  is  strong  enough  to  ignite  the  gas.  In  the 
=11  case  the  transmitted  shock  ignites  the  gas.  However,  the  gas  that  has  passed 
through  the  initial  and  ramp  shocks  is  not  as  highly  compressed  and  heated,  so  it  is  not 
burned  until  it  reaches  the  reflected  shock.  The  transmitted  and  ramp  shocks  ignite  the 
flow  at  Mach  numbers  of  13  and  15;  there  is  a  small  induction  zone  behind  the  ramp  shock 
in  the  =13  case  that  is  not  present  at  Mach  15. 

While  the  methane  mass  fraction  contours  crisply  illuminate  the  burned  and  unburned 
regions  of  the  flow,  they  give  the  impression  that  once  the  fuel  is  burned  there  are  no  more 
chemical  reactions.  Carbon  dioxide  mass  fractions,  shown  in  Figure  5.6,  indicate  this  is 
not  the  case.  After  the  methane  is  consumed,  carbon  atoms  are  distributed  between  CO 
and  C02-  For  hotter  temperatures,  the  balance  shifts  toward  carbon  monoxide;  at 

M„  =  9 

M„  =  11 

M„=13 

M_  =  15 


Figure  5.5  Methane  mass  fraction  contours. 
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relatively  lower  temperatures,  CO^  dominates.  For  the  Mach  7  case,  the  temperature  in 
the  separation  bubble  (behind  the  projectile  shoulder)  is  roughly  32QOK,  and  carbon 
dioxide  concentrations  are  high.  At  increasing  projectile  velocities,  the  temperature  of  the 
gas  in  the  separation  bubble  rises.  For  the  Mach  13  case  the  temperature  is  ~  4000  K,  and 
there  is  practically  no  CO2  •  For  the  higher  speed  cases,  the  state  difference  resulting  from 
the  different  S1/S2  and  S3  shock  strengths  is  visible.  Flow  near  the  tube  wall  is  hotter 
because  it  has  come  through  the  transmitted  shock,  and  there  is  less  carbon  dioxide;  as  it 
expands  downstream,  it  cools,  and  the  COj  fraction  increases. 

The  different  cases  demonstrate  the  ability  of  the  quasi-global  mechanism  to  distin¬ 
guish  subtle  differences  in  gas  conditions,  in  particular  the  detection  of  the  induction  zone 
behind  the  ramp  shock  in  the  Mach  13  case  and  the  state  difference  in  the  Mach  11  case. 
The  mass  fraction  contours  of  carbon  dioxide  and  other  species  (not  presented  here)  estab¬ 
lish  that  all  the  species  included  in  the  gas  model  are  significant  in  some  region  of  the 


M„  =  9 


M„  =  11 


M.  =  13 


M.  =  15 


Figure  5.6  Carbon  dioxide  mass  fraction  contours. 
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flowfield.  This  suggests  that  simulations  including  fewer  species  may  not  be  adequate  for 
ram-accelerator  applications.  In  particular  the  gas  composition  in  high  temperature 
regions  may  suffer. 

As  the  projectile  velocity  increases,  the  ignition  location  moves  upstream.  The  double¬ 
cone  geometry  was  designed  to  delay  combustion  on  the  forebody  by  using  a  shallower 
initial  cone  angle;  the  second  ramp  initiates  combustion,  which  primarily  occurs  on  the 
afterbody.  The  simulations  show  that  through  Mach  15,  at  least,  the  projectile  performs  as 
intended.  There  is  little  combustion  on  the  forebody  upstream  of  the  second  ramp. 
Although  there  is  some  burning  in  the  forebody  boundary  layer,  it  does  not  significantly 
affect  the  flowfield. 

The  peak  temperature  in  the  Mach  15  case  is  approximately  6000  K,  which  is  at  the 
limit  of  the  curve  fits  used  for  calculating  the  equilibrium  constants.  Additionally  the 
assumption  that  nitrogen  does  not  react  is  questionable  at  this  high  temperature.  Nonethe¬ 
less,  simulations  for  =  17  and  =  19  have  been  carried  out.  Mach  number  con¬ 
tours  for  these  flowfields  are  shown  in  Figure  5.7,  and  methane  mass  fraction  contours  are 
presented  in  Figure  5.8.  They  suggest  that  the  initial  shockwave  is  strong  enough  to  ignite 
the  gas  at  Mach  19,  but  not  at  Mach  17.  The  peak  temperatures  in  these  flowfields  are 
observed  along  the  forebody  surface.  Since  nitrogen  is  not  allowed  to  dissociate,  the  tem¬ 
peratures  are  artificially  high.  Consequently  the  simulations  overpredict  the  expansion  of 
the  gas  behind  the  initial  shock.  Despite  this,  unstart  is  not  observed.  If  the  gas  model  is 


M  =17 


M  =19 


Contour  Levels:  Min  =  0.0,  Max  =  19.0,  Increment  =  0.5 


Figure  5.7  Mach  number  contours  for  Mach  17  and  Mach  19  cases. 
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Initial  mass  fraction  of  CH^  =  0.02837 
Contour  levels;  Min  =  0.0,  Max  =  0.028,  increment  =  0.002 

Figure  5.8  Methane  mass  fraction  contours  for  Mach  17  and  Mach  19  cases. 

accurate,  the  simulations  suggest  that  the  projectile  can  be  accelerated  to  Mach  19  without 
disrupting  stable  operation. 

V.4  Discussion  of  the  Separation  Bubble 

Particular  attention  is  now  given  to  the  separation  bubble  generated  by  the  impinge¬ 
ment  of  the  reflected  shock  on  the  projectile  boundary  layer.  For  =  9  case  the  bubble 
is  larger  and  fills  more  of  the  available  tube  cross-section  than  for  the  other  projectile 
speeds;  additionally  the  shape  of  the  recirculation  region  is  not  the  same.  In  this  section, 
the  factors  influencing  the  size  and  shape  of  the  separation  bubble  behind  the  projectile 
shoulder  are  investigated. 

The  separated  region  is  most  easily  seen  in  the  Mach  number  contours  in  Figure  5.3. 
Flow  in  the  bubble  is  primarily  subsonic.  Streamlines  presented  in  Figure  5.9,  for  the 
Mach  11  case,  illuminate  the  region  of  interest.  Several  vortices  can  be  seen  within  the 
separation.  Mach  number  contours  are  also  shown  so  other  flowfield  features  can  be  iden¬ 
tified  and  correlated  to  the  streamline  behavior. 

As  discussed  in  Section  V.3,  the  ignition  location  moves  upstream  as  the  projectile 
velocity  increases.  After  the  gas  is  ignited,  chemical  energy  is  released  into  the  flow.  Con¬ 
sequently  the  temperature  rises.  Since  the  pressure  and  composition  change  only  slightly, 
the  ideal  gas  law  stipulates  the  gas  density  must  decrease.  In  a  flow  system  the  hot  gas 
expands,  filling  a  greater  volume,  and  accelerates.  This  expansion  affects  the  size  of  the 
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Mach  number  contours 


Streamlines 


Figure  5.9  Mach  number  contours  and  streamlines  illustrating  the  separation  bubble. 


separation  bubble. 

For  the  =  7  case,  no  combustion  occurs  except  in  the  forebody  boundary  layer. 
By  the  time  the  burned  gas  reaches  the  recirculation  zone,  any  expansion  due  to  chemical 
energy  release  has  been  assimilated  into  the  flow.  For  projectile  Mach  numbers  of  11,  13 
and  15,  combustion  and  significant  energy  release  take  place  upstream  and  downstream  of 
the  separation  bubble.  However,  for  the  Mach  9  case,  the  reflected  shock  ignites  the  gas. 
Chemical  energy  release  happens  entirely  behind  the  separation.  The  corresponding, 
unbalanced  expansion  enhances  the  shockwave/boundary  layer  interaction  that  generates 
the  separation.  For  the  higher  projectile  speeds,  the  expansion  downstream  of  the  bubble  is 
balanced  by  the  energy  release  upstream. 

To  demonstrate  that  the  chemical  energy  release  intensifies  the  separation  bubble,  the 
Mach  9  simulation  was  repeated  assuming  the  gas  was  inert.  The  composition  is  fixed  to 
the  freestream  proportions.  Figure  5.10  shows  streamlines  for  the  nonreacting  and  reacting 
flowfields.  The  differences  in  the  separation  between  the  Mach  9  case  and  the  higher 
velocities  can  be  seen  by  comparing  Figures  5.9  and  5.10.  For  the  streamlines  starting  near 
the  left  edge  of  each  flowfield,  the  streamline  origins  have  the  same  coordinates.  An  arrow 
marks  the  upstream  limit  of  the  separation  bubble.  The  arrows  and  streamlines  clearly 


Non  reacting 


show  the  separation  is  larger  for  the  reacting  case,  extending  closer  to  both  the  tube  wall 
and  the  projectile  shoulder.  Behind  the  reflected  shock,  the  streamlines  turn  away  from  the 
accelerator  tube  more  severely  in  the  reacting  solution,  indicating  more  expansion. 

V.5  Unstart  Characteristics 

Unstart  was  not  observed  for  any  of  the  projectile  velocities  considered.  The  size  of 
the  separation  bubble  at  Mach  9  is  cause  for  concern,  though.  For  slightly  different 
conditions,  enough  of  the  tube  area  might  be  obstructed  to  choke  the  flow.  In  simulations 
carried  out  by  Soetrisno  et.  al.^  on  canted  projectiles,  a  detonation  wave  propagates 
upstream  when  the  flow  chokes.  The  detonation  wave  passes  the  projectile  tip,  burning  all 
the  fuel  prematurely,  and  the  projectile  accordingly  looses  thrust.  Once  the  detonation 
starts  moving  upstream,  stable  operation  is  lost.  Whether  or  not  the  separation  bubble 
would  unstart  the  flow  is  not  clear,  though.  If  the  flow  choked,  a  detonation  wave  moving 
upstream  would  alter  the  shock  structure  around  the  projectile.  A  change  in  the  shock 
system  might  reduce  the  size  of  the  separation  bubble,  allowing  the  necessary  mass  flow. 
Further  study  is  needed  to  determine  the  role  of  factors  other  than  projectile  speed  on  the 


2.  See  Soetrisno,  Imlay,  and  Roberts  (1993). 


size  of  the  separation,  in  particular  the  gas  composition  and  the  tube  fill  pressure. 

The  double-cone  geometry  was  designed  to  delay  or  eliminate  unstarts  due  to  fore¬ 
body  ignition  of  the  fuel.  From  the  simulations,  it  appears  that  the  geometry  accomplishes 
this  objective.  For  projectile  Mach  numbers  from  7  to  15,  the  initial  shock  did  not  induce 
combustion.  The  simulation  at  Mach  19,  although  questionable  because  of  the  high  tem¬ 
peratures  reached,  suggests  that  even  when  the  initial  shockwave  ignites  the  gas  stable 
operation  is  realized.  Further  modifications  of  the  geometry  or  operational  parameters  (gas 
composition  and  fill  pressure)  may  be  required  to  eliminate  unstarts  caused  by  mechanical 
or  thermal  choking. 


CHAPTER  VI 


CONCLUSIONS 


VI.l  Conclusions 

The  key  findings  of  the  research  presented  in  this  thesis  are  summarized  in  this  sec¬ 
tion.  Remarks  are  made  about  the  chemistry  model  and  its  suitability  for  future  studies. 
Some  statements  concerning  the  modified  diagonal  implicit  scheme  are  made.  Finally,  the 
results  of  the  ram-accelerator  simulations  are  summarized. 

The  most  important  results  of  the  blunt  body  simulations  are  the  findings  concerning 
the  quasi-global  methane  reaction  mechanism.  The  pressure  and  temperature  profiles 
highlight  the  limitations  of  the  model.  In  particular,  the  details  of  the  reaction  front  are  not 
well  predicted.  Temperature  profiles  are  smeared  in  the  energy-release  zone.  As  will  be 
explained  below,  for  ram-accelerator  flowfields  the  smearing  of  the  reaction  front  may  not 
significantly  affect  the  flowfield.  Because  the  single,  global  reaction  represents  a  number 
of  detailed  reactions,  its  range  of  applicability  is  restricted;  however,  with  calibration  to 
established  numerical  or  experimental  data,  the  quasi-global  mechanism  provides  a  rea¬ 
sonable  estimate  of  the  induction  time.  The  main  advantage  of  the  model  is  that  it  is  com¬ 
putationally  less  expensive  than  a  more  detailed  mechanism. 

To  the  author’s  knowledge,  this  thesis  is  the  first  application  of  the  modified  diagonal 
implicit  scheme  to  combusting  flows.  Previous  applications  have  been  for  atmospheric 
reentry  simulations.  A  problem  with  the  scheme  is  that  it  does  not  preclude  the  computa¬ 
tion  of  negative  mass  fractions  for  the  replaced  species.  This  difficulty  has  not  been 
reported  previously.  The  negative  mass  fractions  are  not  the  result  of  discretization  error 
or  floating-point  arithmetic.  Negative  densities  develop  because  the  algorithm  relaxes 
each  species  using  a  different  time  scale,  and  at  the  same  time  elemental  mass  conserva¬ 
tion  is  strictly  enforced.  Even  if  the  negative  mass  fractions  do  not  prevent  the  algorithm 
from  generating  a  solution,  they  may  alter  the  progression  of  the  simulation  in  such  a  way 
that  the  negative  fractions  remain  at  the  final  steady  state.  Different  choices  for  the 
replaced  species  alter  the  evolution  of  the  flowfield,  and  consequently  some  choices  may 
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result  in  negative  mass  fractions,  while  others  do  not.  At  this  point,  careful  testing  through 
trial  and  error  is  the  only  method  of  determining  if  negative  densities  will  arise. 

The  purpose  of  the  research  was  to  investigate  ram-accelerator  flowfields.  Projectile 
Mach  numbers  ranged  from  7  to  15  and  beyond.  As  the  projectile  velocity  was  increased, 
combustion  moved  upstream.  In  most  instances  a  shockwave  separated  the  burned  and 
unburned  gases.  Very  few  induction  zones  were  seen,  and  those  were  very  short.  At  a 
given  speed,  a  shock  that  was  strong  enough  to  ignite  the  gas  generally  did  so  immedi¬ 
ately.  For  this  reason,  the  quasi-global  model’s  smearing  of  the  reaction  front  may  not  be  a 
major  drawback  for  ram-accelerator  simulations. 

The  unstart  characteristics  of  the  double-cone  projectile  geometry  were  sought.  For  the 
conditions  studied,  no  unstarts  were  observed.  For  all  velocities  investigated,  the  gas  in 
the  forebody  boundary  layer  was  burned;  however,  combustion  does  not  spread  to  the  rest 
of  the  flowfield.  Simulations  suggest  that  the  initial  shockwave  ignites  the  flow  above 
Mach  17,  but  that  unstart  does  not  result.  On  the  basis  of  these  two  points,  the  double-cone 
geometry  performs  as  designed. 

An  unstart  appears  most  likely  near  Mach  9.  A  shockwave/boundary  layer  interaction 
results  in  a  large  separated  region  behind  the  projectile  shoulder.  The  size  of  the  separation 
is  greatest  at  Mach  9  where  it  extends  almost  to  the  accelerator  tube  wall.  Unstart  seems 
most  likely  near  this  speed  because,  for  slightly  different  conditions,  the  separation  might 
choke  the  flow.  The  separation  is  largest  at  this  speed  because  combustion  and  energy- 
release  occur  only  on  the  downstream  side  of  the  separation  bubble.  Unbalanced  expan¬ 
sion  of  the  combustion  products  enhances  the  shockwave/boundary  layer  interaction, 
which  enlarges  the  separation. 

VI.2  Future  Work 

Based  on  the  experiences  endured  in  the  research,  several  recommendations  for  future 
investigations  can  be  made.  The  suggestions  are  for  improvements  of  in  two  areas:  physi¬ 
cal  modelling  and  numerical  methods. 

Improvements  in  the  physical  modelling  of  ram-accelerator  flows,  and  hypersonic  and 
combusting  flows  in  general,  can  be  further  divided  into  two  groups.  First,  the  accuracy  of 


the  available  models  can  be  improved.  Specifically,  there  is  no  reliable  chemical  kinetics 
mechanism  for  the  high  pressure  conditions  in  ram-accelerator  flowfields.  A  tremendous 
amount  of  research  has  been,  and  continues  to  be,  devoted  to  hydrocarbon  oxidation; 
unfortunately  the  research  is  concentrated  on  subsonic  flames  at  roughly  atmospheric 
pressures.  Typical  fill  pressures  in  experimental  ram-accelerators  are  from  20  to  50  atm; 
after  passing  through  several  shocks  and  reflections,  pressures  reach  thousands  of  aims. 
The  lack  of  a  reliable  chemistry  model  prevents  accurate  high-pressure  simulations  and 
was  the  major  factor  for  using  1  atm  as  the  freestream  pressure  in  the  simulations  in  this 
thesis.  Perhaps  through  the  use  of  modern  flow  visualization  techniques,  such  as  Planar 
Laser  Induced  Fluorescence  (PLIF),  reliable  kinetics  mechanisms  can  be  developed  for 
high-pressure,  detonation  wave  applications. 

Once  reliable  experimental  data  is  obtained,  numerical  simulations  can  be  used  to  test 
newly  formulated  reaction  mechanisms.  A  technique  for  validating  combustion  chemistry 
models  is  suggested  by  Wilson.^  Blunt-nosed  projectiles  fired  into  detonable  gases  can 
generate  unsteady,  periodic  flowfields,  which  are  highly  sensitive  to  the  kinetics  and  flow 
conditions.^  By  comparing  computational  flowfield  solutions  to  experimentally  measured 
frequencies,  the  performance  of  the  chemistry  model  can  be  assessed. 

Simulations  can  also  be  improved  by  the  inclusion  of  more  accurate  physical  models. 
Despite  the  lack  of  a  validated  chemical  model,  ram-accelerator  simulations  at  higher 
pressures  should  be  carried  out.  Different  fiiel/oxidizer/diluent  combinations  should  also 
be  examined.  While  the  author’s  computational  resources  prohibit  the  inclusion  of  the  pro¬ 
jectile  wake,  current  supercomputers  can  handle  the  additional  requirements.  The  assump¬ 
tions  of  a  thermally  perfect  mixture  and  laminar  flow  could  be  relaxed  with  the  inclusion 
of  the  appropriate  physical  models.  Presently  three-dimensional  and  unsteady  simulations 
are  prohibitively  expensive  when  reacting  gas  models  are  included,  but  with  the  rapid 
advances  in  computer  power  they  are  not  far  in  the  future. 

The  computation  of  negative  mass  fractions  by  the  modified  diagonal  implicit  scheme 


1.  See  Wilson  (1991). 

2.  In  addition  to  Wilson,  numerical  solutions  of  these  unsteady  flows  have  been  obtained  by  Suss- 
man  (1994)  and  Yungster  (1994)  among  others. 


must  be  addressed.  Experience  to  date  suggests  that  the  global  reaction  rate  form  amplifies 
the  difficulties  in  the  current  application.  On  the  other  hand,  the  occurrence  of  negative 
densities  does  not  rely  on  the  inclusion  of  the  global  reaction.  The  conditions  that  lead  to 
negative  mass  fractions  can  be  satisfied  by  a  kinetics  model  comprised  of  detailed  reac¬ 
tions  with  Arrhenius  rate  forms.  Clearly  more  research  must  be  committed  to  this  solution 
advancement  method. 
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THE  ELEMENTAL  FORMULATION 


In  this  section  the  relationship  between  two  sets  of  mass  conservation  equations  is 
made  clear.  Typically  the  species  mass  conservation  equations  are  employed;  one  equation 
is  included  for  each  chemical  species  in  the  gas  model.  This  is  referred  to  as  the  species 
formulation.  In  this  thesis  some  of  the  species  equations  are  replaced  by  elemental  mass 
conservation  equations.^  The  set  comprised  of  ne  elemental  equations  and  ns  -  ne  species 
equations  is  called  the  elemental  formulation. 

Both  formulations  describe  the  same  physics.  Both  implicitly  conserve  total  mass.  The 
species  formulation  implicitly  conserves  elemental  mass,  and  the  elemental  formulation 
implicitly  conserves  the  masses  of  the  replaced  species.  Replaced  species  are  those  whose 
mass  conservation  equation  has  been  replaced  by  an  elemental  conservation  equation.  The 
retained  species  have  conservation  equations  explicitly  included  in  the  set  comprising  the 
elemental  formulation. 

The  conservation  of  mass  is  an  expression  of  the  conservation  of  atoms.  Whether  spe¬ 
cies  masses  or  elemental  masses  are  considered,  the  number  of  atoms  (or  moles)  of  each 
element  must  be  conserved.  The  relationship  between  the  two  formulations  is  based  on  the 
conservation  of  atoms  of  each  element.  The  total  number  of  moles  of  a  given  element  e  in 
the  mixture  is  the  sum  of  the  moles  of  that  element  in  each  chemical  species  s: 


P  e  ^  1  . 

-  —  Cl  1 -  +  ...  Cl  - 


=  z 


e,sM 

s=l  ^ 


(A.1) 


The  species  or  elemental  density  (mass  /  volume)  divided  by  the  molecular  weight  is  the 
number  of  atoms  (moles  /  volume).  The  coefficient  a^  ^  is  the  number  of  atoms  of  element 
e  in  a  molecule  of  species  s.  Consider  carbon  dioxide,  for  example.  For  CO2, 

^c,  CO2  ~  ^  ^o,  CO2  ~  ^  ■ 

The  governing  equations  are  written  in  terms  of  densities.  A  more  convenient  way  of 


1.  The  reasons  for  including  the  elemental  equations  are  discussed  in  Chapter  HI,  as  well  as  in 
Candler  and  Olynick  (1991)  and  Hassan,  Candler,  and  Olynick  (1992). 
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writing  equation  A.l  is  to  include  the  molecular  weights  in  the  coefficients: 


(A.2) 


S  =  I 


where 


^e,  s  ^e, 


(A.3) 


Given  the  constituent  species  of  the  gas  model,  the  coefficients  ^  can  be  determined.  In 
matrix  form,  equation  A.2  can  be  written 


Pi 

«u  •• 

•  «l,n. 

.  \ 

^ne,\  •• 

*  ^ne,  ns 

Pi 

Pns 


(A.4) 


At  times  it  is  necessary  to  determine  the  replaced  densities,  given  the  elemental  and 
retained  densities.  The  coefficient  matrix  can  be  inverted  to  yield 

Pi 


(A.5) 
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P«e+1 

For  the  gas  model  used  in  this  thesis,  there  are  10  species  and  4  elements  {ns  =10  and 
ne  =  4.)  The  elements  and  their  indices  are 

Element:  CHON 

Index:  1234 

The  species  and  their  indices  are 


Species:  CH^  CO^  CO  H^O  OH  O  H 

Index:  123456789  10 


In  this  case  the  first  four  species  are  replaced,  while  the  rest  are  retained.  The  correspond¬ 
ing  coefficient  matrices  are 
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APPENDIX  B 
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PHYSICAL  CONSTANTS 


In  this  appendix  constants  and  parameters  used  in  various  physical  models  are  listed. 
Data  is  from  a  variety  of  sources.  Table  B.l  contains  constants  and  parameters  which 
specify  the  behavior  of  each  species  in  the  gas. 

Collision  integrals  must  be  calculated  to  determine  the  transport  coefficients.  These 
integrals  are  functions  of  the  reduced  temperature,  J* .  For  each  species  the  Lennard- 
Jones  6-12  potentials  are  required  to  determine  the  reduced  temperature.  The  curve  fit 
polynomial  for  the  collision  integrals  is 

f{T*)  =  aQ  +  a{I^  +  a^T*^ +  +  .  (B.l) 


Table  B.2  lists  the  Lennard- Jones  parameters,  while  the  collision  integral  curve  fit  coeffi¬ 
cients  are  given  in  Table  B.3.  Note  that  there  is  a  separate  curve  fit  for  H2O  because  of  its 
polar  nature. 

Curve  fits  are  used  to  specify  the  Gibbs  free  energies  of  each  species  which  are 
required  for  the  determination  of  the  equilibrium  constant  for  each  chemical  reaction.  The 
curve  fit  is  given  by 


(l-lnT) 


C2  -  Co  O  C^  O  Cc  4  Cg 

2  6  12  20  T  i,s 


(B.2) 


Two  sets  of  coefficients  are  listed.  The  first  set  is  valid  for  temperatures  between  300  K 
and  1000  K  and  is  presented  in  Table  B.4.  The  second  set  is  valid  from  1000  K  to  6000  K 
and  is  shown  in  Table  B.5.  Since  nitrogen  is  not  included  in  any  of  the  chemical  reactions, 
its  Gibbs’  free  energy  is  not  required. 

The  chenfical  kinetics  mechanism  is  listed  in  Table  B.6.  Included  are  the  12  reactions 
and  the  associated  parameters  for  Arrhenius  form  reaction  rates. 
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Table  B*l:  Molecular  Constants^ 


Index 

Species,  s 

K 

(kg/kmol) 

7  ° 

K 

(J/kg) 

ndeg(s) 

0 

V,  s,  m 

(K) 

1 _ 1 

1 

CH^ 

16.043 

-4.649153e6 

3.0 

9 

1879.1(3) 

2207.2(2) 

4196.4 

4343.4(3) 

2 

O2 

32.00 

0.0 

1 

2273.7 

3 

H2 

2.016 

0.0 

2.5 

1 

6338.5 

4 

N2 

28.016 

2.5 

1 

3392.1 

5 

CO^ 

44.011 

-8.907252e6 

2.5 

4 

960.14(2) 

1932.2 

3380.3 

6 

CO 

28.011 

-3.930740e6 

2.5 

1 

3121.6 

7 

H20 

18.016 

-1.337154e7 

3.0 

3 

2294.4 

5262.0 

5404.0 

8 

OH 

17.008 

2.286420e6 

2.5 

1 

5374.2 

9 

0 

16.00 

1.551512e7 

1.5 

0 

10 

H 

1.008 

2.154376e8 

1.5 

0 

i.  The  data  is  from  Stull  (1965),  except  for  the  characteristic  temperature  of  vibration  for  OH 
which  is  taken  from  Stull  and  Prophet  (1971). 


Table  B.2:  Molecular  parameters  for  Lennard- Jones  6-12  potentials^ 
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.  Data  is  from  Dixon-Lewis  (1984). 


Table  B.4:  Polynomial  coefficients  for  Gibbs’  free  energies^:  low  temperature  range 
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i.  Data  is  from  McBride,  Heimel,  Ehlers,  and  Gordon  (1963). 

ii.  Valid  for  l(mK<T<6O00K. 


Table  B.6:  Quasi-global  methane  reaction  mechanism^ 


Index,  r 

Reaction 

A“ 

1 

ca  +  0  +  M<^C02  +  M 

5.9xl0'^ 

0 

2055 

2 

0H  +  M<^0  +  H  +  M 

8.0x10'^ 

-1.0 

51,990 

3 

02-^M<^0  +  0  +  M 

S.lxio'^ 

0 

57,650 

4 

7^2  +  M<=>H  +  //  +  M 

2.2x10’'^ 

0 

48,130 

5 

H^O  +  M^OH  +  H-\-M 

2.2x10*^ 

0 

52,640 

6'“ 

C//4+  1/2  O2  027/2  + CO 

4.0x10® 

0 

24,260 

7 

CO  +  OH0CO2  +  7/ 

1.5x10^ 

1.3 

-401.0 

8 

CO  +  O2  <=>  CO2  +  0 

3.1x10" 

0 

18,850 

9 

77  +  02O07/  +  0 

2.2x10’'' 

0 

8422 

10 

1.8x10’® 

1.0 

4462 

11 

2.2xl0’^ 

0 

2557 

12 

H^O  +  O^OH  +  OH 

6.8xl0’^ 

0 

9224 

i.  Data  is  from  Westbrook  and  Dryer  (1981). 

ii.  Units  are  moU  cm,  K,  sec 

iii.  Global  reaction,  original  rate.  In  some  calculations,  is  doubled.  Con¬ 
centration  exponents:  a  =  -0.3  ,  a'  =  -0.5 ,  and  b  =  b'  =  13  . 
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GENERALIZED  COORDINATE  TRANSFORMATION 


The  governing  equations  discussed  in  Chapter  II  are  written  in  index  notation,  which 
is  most  easily  extended  to  Cartesian  coordinates.  However,  a  general  curvilinear  coordi¬ 
nate  system  is  preferred  because  most  interesting  geometries  are  difficult  to  represent  effi¬ 
ciently  in  Cartesian  coordinates.  In  this  appendix  the  transformation  between  Cartesian 
and  general  curvilinear  coordinates  is  delineated. 

The  Cartesian  coordinate  directions  are  x  and  y;  for  the  general  curvilinear  coordi¬ 
nates,  they  are  ^  and  Tj ,  The  chain  rule  of  calculus  is  used  to  relate  derivatives  in  the  two 
systems: 

rrn 

dx  0x3^  3x011  ^  ^ 

and 


0y  0y0^  0y0ll 

where  (p  is  any  property.  The  derivatives  ^  and  ^  are  called  the  metrics  of  the 

transformation.  With  the  abbreviation  for  ^  and  so  on  for  the  other  metrics,  equations 
C.  1  and  C.2  can  be  written  in  matrix  form: 


’0’ 

0X 
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=  [ri 
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0 
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?y. 

0TL 

0TL 

Inversion  yields  an  equation  for  derivatives  in  the  Cartesian  space: 
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_  _ 
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dx 
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dx 
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a^ 

.^y. 

(C.3) 


(C.4) 


where  x^,  x^,  y^,  and  are  referred  to  as  inverse  metrics.  Equation  C.3  can  be 
expressed  in  terms  of  the  inverse  metrics: 
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dx 
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=  [T] 
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3'n. 

(C.5) 


The  inverse  metrics  are  easier  to  compute  than  the  metrics,  so  the  forms  in  equations  C.4 
and  C.5  are  preferred.  The  mesh  metrics  are  used  to  evaluate  the  viscous  terms,  as  dis¬ 
cussed  in  Section  111.4. 

The  surface  cosines  s.' ,  s.' ,  s-' ,  and  s.'  can  be  calculated  by  inspection  of  Figure 

IX  ly  jx  jy 

C.l.  For  the  i  +  1/2  cell  face,  the  surface  is  defined  by 

5,1  = 


(C.6) 


where  the  cell  face  area  is 


>S  2  —  Ax  ^ 


{C.l) 


'■^2 


and  the  surface  cosines  are 


Y 


+  Property  value 
location 

o  Grid  definition 


Figure  C.l  Surface  cosines  for  a  representative  mesh  cell. 
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^y 


fj Ax  ^ 
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-Ax 
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Ax^  +  Ay^ 


Similarly  for  the  j  +  1/2  cell  face, 

^,,1  = 


S  .  =  a/ai^  +  a/ 
^^2 


and 


-Ay 


fj Ax  + 


Ay 


= 


Ax 


Jax^  +  Ay^ 


For  the  i  +  1/2  cell  face, 


A^  = 


Ay  =  y,^i,,.^i-y,.^i,,., 


while  for  the  j  +  1/2  face, 

A^  =  -  hM  Ay  =  y .^1  .^1  -  y . 

The  surface  cosines  and  the  mesh  metrics  are  related  by 


IX 


5. 


and 


= 


^Jy  = 


11, 


ri  2 

a/iIx  +  11, 


(C.8) 


(C.9) 

(C.10) 


(C.ll) 


(C.12) 


(C.13) 


(C.14) 


(C.15) 


Finally,  a  formula  for  the  cell  volume  is  given  by  Hirsch. 

^  =  U (yi,M -yMj)  -  (hj^i (y^j+i-yij) i  ^  (c-i6) 


zind  the  cell  centroid  radius  is  an  average  radius  for  the  cell: 

yc.u  =  i‘-yM.M+yM.i*yi,M*yu'>- 


(.CM) 


1.  See  Hirsch  (1988),  p.  247. 
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APPENDIX  D 


DIAGONALIZATION  OF  THE  FLUX  JACOBIANS 


The  flux  Jacobians  are  required  for  both  the  implicit  and  explicit  fluxes.  The  diagonal 
implicit  scheme  approximates  only  the  inviscid  fluxes,  and  only  the  inviscid  fluxes  are 
considered  in  the  upwind  schemes.  Therefore,  the  fluxes  and  their  Jacobians  discussed  in 
this  appendix  consider  only  the  inviscid  terms.  However,  the  viscous  terms  are  included  in 
the  flowfield  solutions  as  detailed  in  section  III.4,  after  the  discussion  of  the  TVD  method. 

dF'i 

The  flux  Jacobian  in  the  ^  -direction.  A,  is  defined  as  .It  is  calculated  in  this  form 
for  the  implicit  side.  On  the  other  hand,  the  upwind  schemes  (which  are  used  to  evaluate 
the  explicit  side)  require  the  decomposition  of  the  Jacobians.  The  diagonalization  of  A  is 
facilitated  immensely  through  the  use  of  a  similarity  transformation. 

For  convenience,  the  vector  of  conserved  variables,  U,  and  the  vector  of  primitive 
variables,  V,  are  repeated  here: 


U  -  (Pj, ... ,  P„g+p  ••• »  P„^»  pu,  pv,  pE) 

V  =  (Pi,...,p„^,P„^^l,...,p„^,M,V,p)^. 
The  rotated,  inviscid  flux  vectors  are 


Pl«' 

pjV' 

Pne^' 

Pne^y 

p  u' 

r  ns 

G'l  = 

P 

^  ns 

puu'  +  ps-J 

puV+ps.^ 

pvu'+ps.^' 

pvv'+p^.^' 

{pE+p)  u' 

ipE  +  p)  v' 

(D.l) 

(D.2) 


(D.3) 


The  flux  Jacobians  required  for  the  implicit  side  are 


d 


02 


«  o. 

^Q- 


aI; 


w  CX, 

ro  ^ 


to 

ioTl  ^ 


to  ^ 

.  .  .  .  05 

+ 

|q.  1  Q- 

CL 

cS\^ 

rs  ^ 
^  CL 

ro  w 

ro' 

s  ^ 

V 

.K 

& 

co'^ 

- 

j9- 

ro^ 

o 

0 

u 

ro 

V 

cf  ^ 


^  ci,  i 

CL 


s 

^  I 

c  Q- 
ICL  I 


CL  CL, 

C 

cx, 

ro^ 

iCL 

ro^ 

V 

ro 

V 

ro 


idTl 


C  Q. 

iCL  I 


I  .<£?!  Q- 

'  I 


C  CL 

iCL  I 
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d  (pu)  p  d  (pv)  p  L  3  ipE] 


and  5 


ro 

1 

ro 

1 

ro 

V 

ro 

ro 

i>5 

.qTI  Cl. 


«  Q. 
ICL 


+ 

^  Q. 
CL 


50  I 


w  ^  ^ 


04-^ 

50  1 

+ 

n 

50  1 

4 

«  Q. 

CL 

Cl 

iCl  I 

c 

Q- 

CL  1 

ro" 

s  + 


Cl. 

+  Q. 


-  ^ 

\ro  ro 


Iro 

V 

p- 

Cl 

+  CL 


ro  w 


.qTI  c^ 


^  Cl 
c  ^ 

lO, 


p^ 

CD  I 

C  CL 

iCL  I 


P^ 

50  I 

«  Ql 
ICL  \ 


+ 

^  CL 
CL 


r^CL^ 


+ 

+ 

50  Ol 

ro^ 

Cl  ^ 

CL  ro 

ro 

V 

ro 

1 

50 

CL  Cl 

c 

Cl 

ICL 

ro^> 

1 

ro 

4" 


ro 


> 

.qTI  cl 
I 

V 


V 

p»* 

- M 

P 

V 

50  1 

^\  CL 
ICL  \ 

1 

+  ,  .  .  , 

^  CL 

CL 

'•■  4"  ^4- 

1 

^hoT 
^  Iro 

1 

The  derivatives  of  the  pressure  with  respect  to  the  conserved  variables  are 


ne 

i; 


+  (Y-1) 


2  2 
M  +  V 


- 


=  IP, 

r=  1 


R 


■*  .V 


_ 

a(p£) 


In  these  equations, 


Y  =  1  + 


Y-1. 


p7? 


P<^V  +  S"IlP5^V,Vib.. 


(D.6) 


(D.7) 


(D.8) 

(D.9) 

(D.IO) 

(D.ll) 


and 


P^V  ^  Pf^v,  5  ’ 
5=1 


P«=  Ip^' 

5=1 


(D.12) 


(D.13) 


Recall  that  ^  contains  only  contributions  from  translational  and  rotational  energy;  the 
vibrational  energy  contribution  is  handled  separately  in  ^ .  The  definitions  of  ^ 
and  ^  are  given  in  section  II.3  and  11.5,  respectively,  but  are  repeated  here  for  conve¬ 
nience: 


E  =  — 
v,5  M 


ndeg(s) 


^  m=  1 


e. 


V,  s,  m 


re 


expF 


V,  5,  m 


(D.14) 
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"v,  vib,  s 


^v,s  ^  _R 
dT  M 


ndeg  (s)  \ 

I 


r®v,s,/ny_r®v,5. 


expl^- 


^  m  =  1 


exp 


V,  s,  m 


(D.15) 


1 


For  the  explicit  discretization  the  Jacobians  are  decomposed.  The  Jacobian  A  can  be 
diagonalized  to 

-1 


^  =  Pa^Pa 


(D.16) 


where  is  a  diagonal  matrix,  with  the  diagonal  comprised  of  the  eigenvalues  of  A,  and 
the  columns  of  are  the  corresponding  right  eigenvectors.  In  practice  it  is  difficult  to 
find  the  eigenvalues  and  eigenvectors  by  directly  decomposing  A. 

Using  the  identity  matrix  and  the  chain  rule,  the  Jacobian  can  be  expressed  as 


dF'i^  = 

dV  dU  dVdU  dV  dU 


(D.17) 


Defining 


and 


(D.18) 


2  _  dV^F'j 

dU  dV’ 


(D.19) 


equation  D.14  is  rewritten 

A  =  MAM~^ .  (D.20) 

Through  the  use  of  the  similarity  transformation  just  described,  the  determination  of  the 
eigenvalues  is  considerably  simplified.  The  decomposition  of  A  is 

The  eigenvalues  of  A  are  the  same  as  of  A ,  but  the  diagonalization  of  A  is  much  easier  to 
accomplish  because  of  its  simpler  structure.  The  right  eigenvalues  of  A  and  A  are  related 
by 


Pa  =  MP^. 


(D.22) 


In  a  similar  fashion  the  flux  Jacobian  in  the  t]  -direction  can  be  expressed  as 


_ 

dV  dU 


The  matrix  M  is 


dUdVdG'idV  _ 
dVdU  dV  dU  ~ 


(D.23) 


M= 


u  ... 

V  ... 

dpE 

api 

The  terms  in  the  last  row  are 


1  0 
0 

0  ... 


0 

0  1 


u 

V 


1  0  ...  0 
0  ■; 


0 


dpE  dpE 
^Pne  ^Pne+1 


•.  0 
0  1 

poo 

0  p  0 

dpE  dpE  dpE  dpE 
3p  du  dv  dp 


3(p£)  _ 


=  SP 


r=  1 


R  1  j  o  ^  ’ 

re\  \  ^v,r~  j_  j  J  r v,r 


2  2 
M  +  V 


(D.24) 


(D.25) 


3(pg)  _ 

ap 


=  EP 

r=  1 


(D.26) 


du 


ipE)  =  pu 


^(p£)  =  pv 


a(pE) 


1 


(D.27) 

(D.28) 

(D.29) 


dp  y- 1  ’ 

The  matrix  M  ^  can  be  determined  by  inversion  of  M,  but  it  is  simpler  to  calculate  it 
directly.  It  is 


M 


1 

0 


0 


u 

~P 

_v 

~p 

api 


0 


0 

0 

1 


1 

0 


u 
~P 
_v 

~p 

3p«.  ^P; 


dp 


0 


ne+1 


0 


0 

1 

1 

p 

0 

dp  dp 


0 

1 

P 

dp 


0 

0 

dp 


^Pns  ^(P“)  ^(P’')  ^(P^ 


(D.30) 


The  pressure  derivatives  are  the  same  as  those  discussed  previously. 
dF'/  dG'i 

This  matrices  and  are  found  component  by  component.  The  results  are 
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The  matrix  A  is  calculated  by  multiplying  ^  and  : 


A= 


u'  0  ...  0 

0  ■•.  ■•.  : 

:  ■•.  ■•.  0 
0  ...  0  u' 


Pi  Pi " 


ly 


OS.  p 

r  ne  IX  r  ne  ly 


0 


0 


m'  0  ...  0  0 

0  ■•.  ;  •  ;  : 

:  ■•.  ■•.  0  :  :  ; 

0  ...  0  u'  0 

Si' 

u'  0  — 


0 


pa  Si^  pa  Siy  u 


Similarly  for  B , 


B= 


v'  0  ...  0 
0  • 


0 


•.  0 
0  v' 


0 

0 


p,5 


Pi^/v'  0 


1  ^jx  n 


6  s.  '  6  s-  '  0 

^ne  jx  rne  jy 


v'  0  ...  0  P„,.,iVPne+lV  ^ 


'.  0 


Ov^pj'.^  p^.^  0 

r ns  jx  ^ ns  jy 


0 


0 


JX 

p 

^jy 


pa  SjJ  pa  s.y'  v' 


(D.33) 


(D.34) 


The  speed  of  sound,  a,  is  defined  as 


a  =  JyRT ; 


(D.35) 


This  is  not  a  frozen  or  an  equilibrium  speed  of  sound;  the  composition  of  the  gas  is 
accounted  for  in  y  and  R . 

The  eigenvalues  of  A  satisfy 

-^-/J  =  0.  (D.36) 

Solving  this  equation,  the  eigenvalues  are 


detl  A 


f 

u 


/ 

u 

u'+u 

u'-a 


(D.37) 


The  eigenvalue  u'  is  repeated  ns  times;  for  the  equation  set  described  in  the  Chapter  H, 
there  are  ns  +  3  equations  and  thus  ns  +  3  eigenvalues. 

For  the  tj  -direction,  the  result  is 


V  +a 


v'-a 


(D.38) 


The  order  of  the  eigenvalues  is  arbitrary;  however,  the  eigenvectors  should  be  in  the  same 
order  as  the  eigenvalues. 

The  right  eigenvectors  of  A  are  solutions  to 

=  0,  (D.39) 

where  k  is  the  eigenvector  corresponding  to  the  eigenvalue  X.- .  It  turns  out  the  eigenvec- 

A 

tors  of  A  are  not  unique;  there  are  many  arbitrary  constants  allowed.  The  choice  of  these 
constants  can  affect  the  convergence  of  the  numerical  scheme.^  Selections  for  the  arbitrary 
constants  are  made  such  that  the  eigenvectors  are  linearly  independent.  In  this  thesis  the 
following  eigenvectors  are  employed: 


1.  See  Yee(1989). 


r 


and 


2 

a 


2 

a 


Pi 

Pi 

2 

2 

pa 

pa 

1 

^ne 

Pne 

2 

2 

2 

a 

pa 

pa 

1 

Pn^+l  Pne+1 

2 

2 

2 

a 

pa 

pa 

1 

Pn. 

Pn. 

2 

2 

2 

a 

pa 

pa 

o  ' 

_r  ' 

Jy. 

^ix 

pa 

pa 

pa 

V 

pa 

pa 

pa 

0 

1 

1 

Pi 

Pi 

2 

pa 

2 

pa 

1 

Pne 

Pne 

2 

a 

2 

pa 

2 

pa 

1 

P«e+1  Pne+l| 

2 

a 

2 

pa 

2 

pa 

1 

Pn5 

P  ns 

2 

a 

2 

pa 

2 

pa 

/ 

V 

pa 

pa 

pa 

~^Jx 

V 

pa 

pa 

pa 

0 

1 

1 

(D.40) 


(D.41) 
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These  matrices  are  inverted  to  obtain  P-  and  P- 

A  B 


and 


(D.42) 


(D.43) 


The  eigenvectors  of  A  and  B  are 


Pa  = 


a 


a 

V 


a 

l^(pE) 
^api 


a 


Pb  = 


a 


a 

IdipE) 


^Pi 


fL 

K  ■ 

pa^ 

pa^ 

Jl_ 

Pne 

Pne 

pa^ 

pa^ 

1 

Pne+1 

Pne+1 

2 

^  2 

2 

a 

pa 

1 

P/w 

P  ns 

pa^ 

pa^ 

u 

V 

u-as.; 

u  +  a5j.^' 

a 

V 

V 

1 

v  +  as./ 

a 

2 

a 

IdipE) 

1  a(p^ 

ia(pE)  -v' 

pE^p  u' 

pE+p  1  m' 

^Pne 

a2^P„,+i  ■■ 

■  a^^Pns  « 

pa^  ^ 

pa^  « 

h. 

A 

2  2 

pa  pa 


J_ 

1 

1 

Pne 

2 

pa 

pa^ 

1 

Pne+l 

Pne+l 

a^ 

pa^ 

pa^ 

1 

a2 

pa^ 

2 

u 

V 

M  -  as.' 

M  +  ai'.  ' 

JX 

a 

V 

-V 

v-as.; 

V  +  asjy 

a 

ia(p^ 

1  a(p^ 

ia(p^ 

u' 

pE+p  v' 

pE+p  ^  v' 

a^^Pne+l 

■  a^3p„ 

a 

pa^  ^ 

pa^  ^ 

(D.44) 


.(D.45) 
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